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CHAPTER  I 


INTRODUCTION 


This  is  the  second  complete  documentation  of  the  HELP^^^ 
code  since  it  was  developed  five  years  ago.  For  those  un- 
familiar with  HELP,  the  following  remarks  will  give  some 
historical  perspective.  The  HELP  code  evolved  from  three 
major  hydrodynamic  codes  developed  over  the  past  twenty  years: 
RPM,  OIL,  and  PIC.  Starting  from  the  one-material  code, 
RPM^^^,  the  HELP  code  added  the  capability  of  describing 
multimaterial  interactions  and  of  treating  material  strength 
as  elastic-plastic  rather  than  rigid-plastic.  RPM,  in  turn, 
was  an  extension  of  OIL^^^,  a compressible  fluid  code,  and 
included  material  strength  in  the  rigid-perfectly  plastic 
approximation.  The  OIL  code  was  derived  from  the  PIC^^^  code 
and  replaced  the  discrete  PIC  particles  by  a continuum.  This 
latter  transition  to  a continuum  calculation  was  made  in  order 
to  permit  the  treatment  of  very  small  compressions  without 
having  to  use  an  extremely  large  number  of  particles  for  ade- 
quate mass  resolution.  The  disadvantage  of  this  transition, 
however,  was  the  loss  of  certain  Lagrangian- type  features  of 
PIC,  and  specifically  the  ability  to  treat  flows  containing 
more  than  one  mate'll  al.  The  m;  ;or  purpose  of  developing  the 
HELP  code  was  to  restore  che  multimaterial  capability  of  PiC 
while  retaining  the  ease  of  treating  small  compressions  and 
the  computing  economy  of  a continuous  Eulerian  model. 

The  accurate  treatment  of  material  interface  or  free 
surface  motion,  however,  requires  a detailed  knowledge  of 
surface  location, which  is  not  contained  in  a purely  Eulerian 
model.  The  developers  of  the  HELP  code,  therefore,  deviated 
from  the  pure  Eulerian  model  by  employing  massless  tracer 
particles  to  track  material  interfaces  and  free  surfaces. 
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These  particles  are  initially  placed  along  the  surface  bound- 
aries and  are  moved  thereafter  with  a local  material  velocity, 
thus  creating  a Lagrangian-type  definition  of  moving  surfaces 
without  sacrificing  the  capability  of  easily  treating  extreme 
distortions.  Furthermore,  diffusion  of  materials  across 
these  surfaces  is  prevented  in  HELP  by  taking  into  account 
the  surface  positions  when  transporting  material  across  the 
Eulerian  cell  boundaries. 

In  the  years  since  HELP  was  first  developed  and  docu- 
mented it  has  been  applied  to  a wide  variety  of  multimaterial 
problems  in  continuum  mechanics.  In  the  process  it  has  been 
significantly  expanded  to  increase  its  generality  and  applic- 
ability. For  example,  sliplines  were  added  in  order  to  model 
more  accurately  the  plugging  failure  in  thin  targets^^’^^ 

r 7 1 

and  the  high-explosive  metal  interactions  in  shaped  charges  ^ 
and  fragmenting  munitions . Furthermore,  a high  explosive 
detonation  model  was  added  which  would  account  for  multiple 
explosives,  multiple  detonation  points  and  wave  shapers. 

As  the  applications  of  the  HELP  code  increased  the  need 
for  a more  general  problem  generator  grew.  The  original 
problem  generator  was  written  primarily  for  calculations 
involving  spheres  and  cylinders  impacting  finite  and 
semi-infinite  plates.  This  original  generator  clearly  was 
inadequate  for  generating  shaped  charge  problems  or  even 
projectiles  with  various  shaped  noses.  The  new  HELP  gener- 
ator can,  without  modification,  generate  any  configuration 
where  the  package  boundaries  can  be  decomposed  into  straight 
line  segments  or  arcs  of  circles  and  ellipses. 

In  addition  to  these  improvements,  some  aspects  of  the 
code  have  been  reformulated.  For  example,  the  hydrodynamic 
phase  (HPHASE)  has  been  simplified  (the  equations  are  solved 
in  a single  pass)  in  order  to  provide  more  accurate 
definitions  of  pressures  and  velocities  near  free  surfaces. 
(See  Section  2.2.2.)  This  has  eliminated  the  need 
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to  "glue"  free  surface  cells,  a method  v/hich  proved  unsatis- 
factory in  certain  situations.  The  single  pass  formulation 
also  simplified  the  addition  of  an  artificial  viscosity  term 
to  cell  boundary  pressures.  Another  change  involved  the 
ordering  of  the  three  phases.  It  has  been  found  that  more 
realistic  strain  rates  are  computed  using  the  post  - transport 
cell  velocities.  Therefore  the  strength  phase  now  precedes 
the  hydrodynamic  phase  in  the  computational  cycle.  This 
reordering  of  the  phases  also  led  to  their  renaming:  PH3  has 

become  SPHASE;  PHI  has  become  HPIIASE,  and  PH2  has  become 
TPHASE.  And  finally,  the  INFACE  subroutine  has  been  broken 
down  into  five  subroutines,  thereby  reducing  the  core  memory 
requirements  of  the  code  when  it  is  segmented,  and,  perhaps 
more  importantly,  clarifying  the  use  and  treatment  of  the 
material  interfaces. 

A summary  of  the  major  chapters  of  the  text  follows. 
Chapter  II  gives  a general  description  of  the  numerical 
method  and  the  material  models  employed  by  HELP.  The  grid 
boundary  conditions  are  discussed  in  detail  in  Chapter  III 
as  they  apply  to  each  phase  of  the  calculation.  Chapter  IV 
offers  a discussion  of  the  material  interfaces,  as  well  as 
the  determination  of  the  pressure  and  shear  yield  strength 
for  multimaterial  cells.  The  equations  which  govern  the 
slipline  cells  are  given  in  Chapter  V along  with  general 
guidelines  on  their  use,  and  in  Chapter  VI  the  mechanisms 
added  to  the  code  for  modelling  plugging  failure  are  described. 
Procedures  for  generating  and  restarting  a HELP  calculation 
are  given  in  Chapter  VII  as  well  as  some  general  statements 
about  the  code's  capabilities  and  limitations.  The  inter- 
active options  available  to  the  user  are  described  in  Chapter 
VIII,  and  a list  of  the  error  conditions  are  discussed  in 
Chapter  IX.  In  Chapter  X a general  flow  diagram  of  the  code, 
a description  of  all  the  subroutines,  and  a dictionary  of 
variables  is  given.  Several  sample  input  decks  for  the 
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various  types  of  problems  to  which  HELP  is  commonly  applied 
are  given  in  Chapter  XI.  The  HELP  output  is  also  described 
in  this  chapter.  The  last  chapter,  Xll,  describes  several 
applications  of  the  HELP  code.  The  three  appendices  give 
instructions  for  dimensioning  the  arrays,  an  input  form 
with  abbreviated  instructions  for  generating  initial  conditions, 
and  a diagram  indicating  how  the  code  might  be  segmented  to 
save  core  storage. 

This  second  documentation  of  the  HELP  code  is  consider- 
ably more  detailed  than  the  original  and  is  more  oriented 
toward  the  new  user  of  the  code.  As  with  any  large,  complex 
computer  code,  HELP  is  most  effectively  used  by  those  persons 
who  have  a good  grasp  of  its  numerical  methods,  material 
models,  conventions,  assumptions,  capabilities  and  limitations. 
It  is  hoped  by  the  present  authors  that  this  document  will 
provide  this  basic  information.  Any  suggestions  for  improve- 
ments in  the  code  or  its  documentation  will  be  appreciated 
by  the  authors. 
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CHAPTER  II 

GENERAL  DESCRIPTION  OF  THE  NUMERICAL  METHOD 


In  the  present  section  the  conservation  equations  and 
the  finite  difference  analog  to  treat  these  equations 
are  given.  The  material  model,  which  includes  strength 
effects,  is  also  presented. 


I 


2.1  CONSERVATION  EQUATIONS 

Space  is  divided  into  fixed  Eulerian  cells  through 
which  the  material  moves.  To  arrive  at  expressions  for  the 
rate  of  change  of  total  mass,  momentum  and  energy  within 
such  a cell,  it  is  convenient  to  start  with  the  conservation 
equations  governing  the  motions  and  interactions  of  contin- 
uous media  in  the  form: 


f?-'  ( 

H 

(Continuity  Equation) 

^ (pu.) 
3t  dx7 

(2.1) 

1 

Du.  - 

! 

f 

(Equation  of  Motion) 

P Tt  " -577 

(2.2) 

II 

r > 

(Energy  Equation) 

DEt  g 

^ UT  " ^ • 

(2.3) 

Here  is  the  stress  tensor,  which  can  be  regarded  as  the 
sum  of  the  hydrostatic  stress,  and  a stress  deviator 
tensor,  s^j , i . e. , 


. . = s. . - 6. .P 
11  11  11 


and  ^i^i  ^i  total  energy  (kinetic  plus 

internal)  per  unit  mass.  Tensor  notation  is  implied,  so 
that  repeated  indices  denote  summations. 
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Expanding  the  convective  derivatives  in  Eqs.  (2.2)  and 
(2.3),  Df/Dt  = 9f/3t  + u.  3f/3x.,  then  adding  Eq.  (2.1)  times 
u.  to  Eq.  (2.2),  and  F.q.  (2.1)  times  to  Eq.  (2.3),  and 
collecting  terms,  gives 


^ ^ “ij  ■ -sIt 


(2.5) 


(2.6) 


For  the  developments  to  follow  it  is  desirable  to  re- 
place these  differential  equations  by  the  analogous  integral 
equations,  obtained  by  integrating  over  the  cell  volume,  V, 
and  then  converting  the  volume  integral  of  divergences  to 
surface  integrals  over  the  cell  surfaces.  Equations  (2.1), 
(2.5),  and  (2.6)  then  become 


Jt 


f pdv  - - / pUj^n^dS 


(2.7) 


3 

Jt 


f pu.dV  = J Oijn.dS  J puju.njds  (2.8) 


3 


f pE.j,dV  =J  o^jU^n^dS  -J  pu^E.j.n^dS  . (2.9) 
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II 


2.2  DIVISION  INTO  PHASLS 

It  is  convenient  to  express  the  integral  conservation 
relations,  Eqs.  (2.7)  through  (2.9),  as  finite  difference 
equations  over  the  time  step  At  and  also  to  decompose  the 
total  stress,  o^j,  into  its  deviator  and  hydrostatic  components, 
according  to  Eq . (2.4).  This  gives,  for  the  increments  of 
total  mass  (m) , momenta  (muj)  and  energy  (mE.p)  within  the 
cell, 


Am  = 


A(mUj) 


A(mE.p) 


tf  Sj.njdS 

■'s 

tf  Sj.u.njdS 


PnjdS 


Pu^n^dS 


t / pu.n.dS 

•'s  (2 

tf  (pu.u.)nj 

S 

(2 

tj  (pu^E.p)n^ 


(2.10) 


(2.11) 


(2.12) 


Here,  the  terms  on  the  right  are  divided  into  increments  due 
to  stress  deviator  forces  on  the  cell  surface  (first  column), 
those  due  to  tlie  pressure  forces  on  the  cell  surface  (second 
column),  and  those  due  to  the  transport  of  mass,  momentum  and 
energy  through  the  surface  of  the  cell  (third  column).  These 
three  types  of  increments  are  accounted  for  in  distinct  phases 
of  the  computation.  Specifically,  during  each  time  step  all 
cells  are  updated  three  times  to  account  for  the  following: 

• Effects  of  the  stress  deviators  (SPHASE,  Strength  Phase)^ 

• Effects  of  pressure  (HPHASE,  Hydrodynamic  Phase)^ 

•Effects  of  transport  (TPHASE,  Transport  Phase). 

The  phases  are  calculated  in  the  order  listed.  In  the  dis- 
cussion that  follows,  the  calculation  of  the  terms  on  the 
right  of  Eqs.  (2.10)  through  (2.12)  are  described  sequentially 
starting  with  SPHASE.  The  equations  are  written  in  cylindrical 

* See  DT  in  Section  10.3. 


coordinates  although  tlie  code  models  both  the  axisynunetric 
and  plane  strain  cases. 

Some  preliminary  definitions  will  be  useful.  Super- 
script n on  a variable  refers  to  the  value  of  the  variable 
at  the  beginning  of  the  time  step  and  superscript  (n+1)  de- 
notes the  value  at  the  end.  In  this  discussion  a typical 
cell  in  the  interior  of  the  grid  is  considered;  Chapter  III 
has  a discussion  of  the  special  conditions  which  describe 
the  calculation  at  the  grid  boundaries  or  at  the  axis  of 
symmetry.  For  a typical  cell,  denoted  by  a value  of  the 
index  k,  the  dependent  variables  for  that  cell  are  written 
P(k),  u(k),  v(k) , Ej(k),  m(k),  representing  respectively 
the  pressure,  radial  and  axial  components  of  velocity,  the 
specific  internal  energy  and  the  mass  for  cell  k.  The 
adjacent  cells  above,  below,  to  the  right  and  left  of  k 
will  be  designated  respectively  as  ka,  kb,  kr  and  kn . Here 
the  terms  above,  below,  right  and  left  refer  to  a cross- 
section  view  of  the  cells,  in  which  the  left  border  is 
parallel  to  the  axis  of  symmetry  with  z increasing  upward 
(see  Figure  2.1).  Each  cell  is,  then,  the  torus  obtained 
by  rotating  the  rectangle  (since  Ar  f Az  in  general)  about 
the  axis  of  symmetry. 

2.2.1  SPHASE  - The  Effects  of  Deviatoric  Stresses 

This  section  discusses  the  modification  of  the  cell 
velocity  and  internal  energy  due  to  the  stress  deviators. 
The  derivation  of  the  stress  deviators  is  given  in  Section 
2.3.2.  The  SPHASE  subroutine  should  be  executed  only  if  at 
least  one  material  in  the  problem  has  shear  yield  strength. 

2. 2. 1.1  Continuity  Equation.  (2.10) 


No  contribution  in  SPHASE. 


2. 2. 1.2  Equation  of  Motion,  (2.11' 


The  SPHASn  contribution  to  the  equation  of  motion  (2.11) 


f 


Ag(muj)  = At 


f s.  n 

7 1)  a 


dS  - 


The  axial  and  radial  momentum  increments  will  be  discussed 
separately. 


Effect  on  Axial  Motion  - In  this  case  the  velocity  Uj  is 
the  axial  velocity,  and  the  term  s^^n^  is  the  axial  component 
of  the  deviatoric  stress  on  a surface.  For  a torus  with  rectang- 
ular section,  s^jn^  on  the  various  faces  is 


s^^^  cell  top 


- s^^^  cell  bottom 


s^,-'  outer  cell  surface 
rz 


- inner  cell  surface 


where  zz  and  rz  subscripts  denote  normal  and  shear  stresses. 
The  top,  bottom,  right  and  left  surfaces  are  indicated  by 
(t) » (b) , (r),  and  ()i),  respectively.  The  equation  for  axial 
motion  is  therefore 


Ag(mv) 


V( 


t\  - r^lAt  + s^^^  27rr2  Az  At 


- s^*'^  27rri  Az  At 


where  rj^  and  r2  are  the  inner  and  outer  radii,  respectively, 
of  the  cell  and  Az  is  its  axial  dimension. 


..'ibV''"  -'-i; 


‘ ' i ■'  I'iiiiftfttolft 
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In  this  equation,  and  in  the  remainder  of  this  section, 
cell  boundary  stresses  are  determined  from  a simple  average 
of  the  time  n cell  centered  stresses,  e.g.. 


(t)  . 


(r)  . 


Effect  on  Radial  Motion 

In  this  case  the  velocity  u^  is  the  radial  velocity,  and 
the  radial  components  Sj^jn^  of  the  deviatoric  stresses  on  the 
various  faces  of  the  mass  element  as  shown  in  Figure  2.2  are 


s^  •'  top  of  mass  element 
zr  ‘ 


bottom 

zr 


right 


left 

rr 


- Sgg(A0/2)  front 

- Sfl„(A0/2)  back  . 


The  mass  of  the  element  is  mA0/2Tr,  where  m is  the  cell  mass. 
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17  ^s 


A,  (mu) 


4V)  (4  - 4 


^ At  + T2  ^6  Az  At 


- r,A0  Az  At  - Sq^AzAt  A6  At 

rr  1 06 

or  (multiplying  by  27t/A0)  the  equation  for  radial  motion  is 


Ag(mu) 


= 2tt  At 


H - ^i) 


♦ s(''> 


s;_-'  r-  Az 
rr  L 


rj^  Az  - SqqAz  ArJ 


2 . 2 . 1 , 3 Energy  Equation,  (2.12) 

In  the  SPHASE  term  from  Eq.  (2.12), 

Ag(mE.j,)  = Aty  s..u^n.dS 

s 

Si^UjUi  is  the  work  rate  per  unit  surface  area  which,  for 
the  various  faces  of  the  torus,  is 

(s(^)v^^)  + cell  top 

\ Z Z Z IT  f 

cell  bottom 

\ zz  zr  / 

/s(r)y^^^  + outer  cell  surface 

y rz  rr  / 

- inner  cell  surface 

\ rz  rr  / 

where  (t) , (b) , (r) , (i)  denotes  stresses  and  velocities  at 
the  top,  bottom,  right  and  left  cell  faces,  respectively.  The 


cell  face  velocities  are  determined  from  a simple  average  of 
the  time  n velocities  in  the  two  cells,  e.g., 


(t)  , v"(k)  ^ v"(ka) 
2 


, u"[k)  ^ u"(klt) 
2 


The  equation  for  the  SPHASE  change  in  total  cell  energy  is 


Ag(mE, 


3(t)^(t) 

zz 


zr 


(b)  . 


3(b)uCb: 

zr 


m - ^i) 


Ws(r)^(r)  , 3(r)„(r)\  3,  4^  4^ 

\ rz  rr  / z 

■ ('rz^'^^^^  * ^rr^^^^^)  ’^l 

The  post-SPHASE  cell  specific  kinetic  energy,  1/2  u^u^, 
is  determined  by  the  post-SPHASE  velocities  obtained  from 
the  equation  of  motion.  The  post-SPHASE  specific  internal 
energy  for  the  cell,  Ej,  is  then  calculated  from  the  updated 
value  of  specific  total  energy  and  specific  kinetic  energy. 


Ej  = E.J.  - 1/2  u^u^ 

The  partitioning  of  the  internal  energy  increment  among  the 
various  materials  in  a multimaterial  cell  is  discussed  in 
Section  4.6.3. 

This  completes  the  discussion  of  the  momentum  and  energy 
increments  due  to  the  deviator  stresses  as  calculated  in 
SPHASE. 
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2.2.2  HPHASE  - The  Effects  of  Pressure 

This  section  discusses  the  modification  of  the  cell 
velocity  and  internal  energy  due  to  hydrostatic  pressures. 
The  discussion  will  be  in  terms  of  pressures  acting  on  cell 
boundaries.  Sections  2. 2. 2. 4 and  2.2. 2.5  will  discuss  the 
definitions  of  these  cell  boundary  pressures. 


2. 2. 2.1  Continuity  Equation,  (2. 10] 
No  contribution  in  HPHASE. 


2. 2. 2. 2 Equation  of  Motion,  (2.11) 

The  HPHASE  contribution  to  the  equation  of  motion  (2.11) 


A„(mu.)  = -At 


/ ""j' 


The  axial  and  radial  n.omentum  increments  will  be  discussed 
separately. 


Effect  on  Axial  Motion 

In  this  case  u^  is  the  axial  velocity  and  nj  = -1,  0,  +1 
are  the  axial  components  of  the  unit  normal  to  the  bottom, 
sides  and  top  of  the  cell,  respectively.  The  equation  of 
axial  motion  is  therefore 


A|_j(mv)  “ P^^ff|r2  - r^jAt  - p'’'*V|: 


r‘  - rjAt 


where  r.,  and  r,  are  the  radii  of  the  outer  and  inner  cell 
surfaces,  respectively,  and  P^  , P^  ^are  the  pressures  on  the 
bottom  and  top  cell  faces. 
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Effect  on  Radial  Motion 

In  this  case  the  velocity  Uj  is  the  radial  velocity. 
Referring  to  Figure  2.2,  the  radial  components  of  the  unit 
normal  take  the  values  n^  = -1,  1,  -A6/2,  -A0/2,  0,  0, 
corresponding  to  the  inner  and  outer  surfaces,  the  two  side 
surfaces  and  the  top  and  bottom,  respectively.  The  pressure- 
integral  contribution  in  Equation  (2.11)  becomes 

^ Aj^(mu)  = P^^r^AS  Az  At  - P^*'1'2A6  Az  At  I^^]\r  Az  A0  At 

where  P^^^^and  P^^^re  the  left,  right  and  side  face 

(s) 

pressures,  respectively.  The  side  pressure  P is  taken  to 
be  the  average  of  the  left  and  right  boundary  pressures 


f!^>. 


Since  Ar  * r2  - r^  the  equation  for  radial  motion  becomes 


Ap^(mu)  = '2  ZiTAzAt  . 


2. 2. 2. 3 Energy  Equation,  (2.12) 

In  the  HPHASE  term  from  Equation  (2.12), 

Aj^(mE^)  = -At  J Pu^n^dS  , 

S 

Pu^h^  is  the  work  rate  per  unit  surface  area  which,  for  the 
various  faces  of  the  torus,  is 
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p(a)^(a) 

.pCb)^(b) 

pCr)u(r) 

.p(0^(0 


cell  top 
cell  bottom 
outer  cell  surface 
inner  cell  surface 


where  the  superscripts  a,  b,  r,  t.  denote  pressures  and  veloc- 
ities at  the  top,  bottom,  right  and  left  cell  faces,  respectively. 
Multiplying  these  terms  by  the  time  step  and  by  the  surface 
area  they  are  acting  upon  gives  the  HPIIASE  equation  for  the 
change  in  cell  total  energy  as 


r I 


A^(mE^)  = - AtKP^^^v^^^  - 7T  (r^  - r^)  + 

(p(r)u(r)r^  - P^*"^u^^^rj)  2ttAz]  . 

The  post-HPHASE  cell  specific  kinetic  energy,  1/2  u^u^,  is 
determined  by  the  post-HPHASE  velocities  obtained  from  the 
equation  of  motion.  The  post-HPHASE  specific  internal  energy, 
Ej,  is  then  calculated  from  the  updated  value  of  specific  total 
energy  and  specific  kinetic  energy: 


F = F 


u.Ui 


The  partitioning  of  the  updated  internal  energy  among  materials 
in  multimaterial  cells  is  discussed  in  Section  4,6.3. 

2 . 2 . 2 . 4 Treatment  of  Free  Surfaces  and  Large 
bensity  biscbntinuities 

Thus  far,  the  HPHASE  equations  for  momentum  and  energy 
increments  have  been  derived  in  terms  of  cell  boundary  pres- 
sures and  velocities.  Normally,  the  cell  boundary  values 
are  defined  to  be  the  simple  average  of  the  post-SPHASE  cell 
centered  quantities  on  either  side  of  the  boundary,  e.g., 


P(k)  P(kb) 


P(k)  + P(ka) 
,(r)^  u(k)  ^ u(kr) 


,U)^  u(k)  •»•  u(k&) 


However,  these  definitions  can  lead  to  excessive  acceleration 
of  partially-filled,  free-surface  cells  or  of  cells  containing 
relatively  low  density  material.  In  the  case  of  a full  cell 
adjacent  to  a partially- filled,  free- surface  cell,  it  would 
seem  that  the  condition  at  the  boundary  between  them  should 
relate  more  closely  to  the  void  than  the  normal  averaging 
procedure  would  allow.  The  over-acceleration  of  the  partially 
filled  cells  would  be  prevented  by  reducing  the  pressure 
gradient  acting  on  them. 

One  method  of  requiring  a boundary  pressure  to  reflect 
the  presence  of  voids  and  low-density  materials  more  closely 
is  to  define  the  cell  boundary  pressure  to  be  an  inverse 
compression  weighted  average  of  the  adjacent  cells,  i.e., 
for  cells  k and  kb  the  boundary  pressure  between  them  is 
defined  to  be 

C(kb)P(k)  + CCk)P(kb) 

^ c(ki)T'+  rcki — 

where  C(k)  and  C(kb)  are  defined  (for  pure  cells  and  free- 
surface  cells  having  only  one  material)  as  follows: 

C(k)  = [m(k)/VOL(k)]/p^ 

C(kb)  = [m(kb)/VOL(kb)]/p„ 


2-14 


whore  m(k)  , m(kb) , VOL(k),  VOL(kb),  are  the  total  masses  and 
total  volumes  of  their  respective  cells,  and  p is  the  reference 
density.  The  above  compression  definitions  are  meant  to  be 
used  as  weighting  factors,  and,  in  the  case  of  one  material 
free-surface  cells,  assume  the  material  is  smeared  over  the 
entire  volume  of  the  cell.  For  multimaterial  cells  the  com- 
pression term  is  a simple  average  of  the  compressions  of  the 
materials  in  the  cell. 

The  effect  of  inverse  compression  weighting  is  to  weigh 
more  heavily  the  pressure  associated  with  the  underdense  cell 
when  defining  the  cell  boundary  pressure.  This  lowers  the 
pressure  gradient  acting  on  the  underdense  cell.  It  should 
be  noted  that  if  the  neighboring  cells  have  the  same  compression, 
inverse  compression  weighting  reduces  to  the  simple  averaging 
procedure  used  before. 

For  boundaries  where  inverse  compression  weighting  is 
used  to  determine  the  cell  boundary  pressures,  the  cell 
boundary  velocities  used  in  the  energy  equation 

are  defined  to  be  a compression  weighted  average,  e.g., 

„(b)_  cmufk)  + £(kb)u(kb) 


„(b)_  C(k)v(k)  + C(kb)vCkb) 

" C(k)  +”C(kb)  ' • 

This  definition  weighs  more  heavily  the  velocity  of  the  denser 
cell  and  assumes  its  velocity  is  probably  more  realistic  than 
the  velocity  of  the  less  dense  cell. 

The  process  of  inverse  compression  weighting  for  boundary 
pressures  and  compression  weighting  for  boundary  velocities 
is  used  at  all  boundaries  of  a free-surface  cell. 

However,  since  the  over-acceleration  problem  can  occur  when- 
ever a dense,  high-pressure  cell  is  adjacent  to  an  underdense. 
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low-pressure  cell,  an  input  variable,  CRATIO,  has  been  added 
to  the  code  in  order  to  allow  the  user  to  control  the  use  of 
the  compression  weighting.  Whenever  the  compression  ratio 
between  two  cells  is  greater  than  CRATIO,  the  pressure  at  the 
boundary  between  them  is  taken  to  be  the  inverse  compression 
weighted  average  of  the  cell-centered  values,  and  the  boundary 
velocity  is  taken  to  be  the  compression  weighted  average. 

4 

CRATIO  has  a default  value  of  10  which  essentially  restricts 
the  use  of  the  weighted  averages  to  the  free  surface. 


2.2.2. 5 Artificial  Viscosity 

The  stability  of  a HELP  calculation  depends,  in  part, 
on  the  effective  viscosity  which  results  from  converting 
kinetic  energy  into  internal  energy  during  mass  transport. 

If  velocities  are  small,  there  is  little  mass  flux  between 
cells,  and  as  a result,  little  effective  viscosity.  Since 
these  conditions  have  arisen  in  a variety  of  calculations, 
an  option,  LVISC,  has  been  added  which,  when  set  equal  to 
one,  will  signal  the  code  to  add  an  artificial  viscosity 
term  to  the  cell  boundary  pressures  in  subroutine  HPHASE. 

The  discussion  which  follows  presents  a form  of  artificial 
viscosity  which  has  been  useful  in  several  calculations.  How- 
ever, when  the  code  is  applied  to  a new  class  of  problems,  it 
may  be  that  other  artificial  viscosity  formulations  will  give 
improved  results.  The  user  can  make  such  changes  if  the  need 
arises . 

fa) 

The  definition  of  the  artificial  viscosity,  Q'  , which 
can  be  added  to  the  pressure  at  the  top  cell  boundary  is 


( v(k)  - v(ka)).  ^ 


[ 

P (k)  + P (ka) 

. P(k)  V.  P(ka) 

i 

2' 

2 
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where  p is  material  density. 

An  analogous  term,  added  to  the  pressure  at  the 
bottom  cell  boundary,  can  be  defined  by  substituting  kb  for 
ka  in  the  above  equation. 

For  the  right  cell  boundary,  the  artificial  viscosity 


* 1 11  (!<]  ■ 11  (kr)j  . 


P fkl  + P (kr; 


p(k1  + P(kr: 


(j,) 

and  an  analogous  term,  Q , is  defined  for  the  left  cell  boundary 
by  substituting  kil  for  kr. 

If  one  of  the  cell  boundaries  is  also  a transmittive  grid 
boundary,  the  artificial  viscosity  term  for  that  cell  boundary 
is  zero.  If,  instead,  one  of  the  cell  boundaries  is  a reflec- 
tive grid  boundary,  the  artificial  viscosity  for  that  cell 
boundary  is 

-2  . o)(k)  Vl  P(JOl  • P (kT  , 
where  co  is  the  velocity  component  normal  to  that  boundary. 


2.2.3  TPHASE  - The  Effects  of  Transport 


The  purpose  of  this  section  is  to  describe  how  the  trans- 
port of  mass,  momentum  and  energy  from  cell  to  cell  is  accounted 
for  in  the  code.  This  is  done  by  calculating  the  integrals 
in  the  last  terms  of  Eqs.  (2.10)  through  (2.12).  In  the  dis- 
cussion below,  attention  is  given  to  the  case  of  transport 
between  pure  cells.  The  special  provisions  required  to 
treat  interface  cells  are  described  in  Chapter  IV. 
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2*2. 3.1  Continuity  Equation.  C2.10] 

From  Eq.  (2.10)  the  mass  being  transported  is 


A^(m)  = - At  y*  pu^n^dS 


and  is  determined  for  each  cell  face  from 


6m  = - u^A^At 


where  p°  is  the  density  of  the  cell  from  which  the  mass  moves 
(donor  cell),  is  the  area  of  the  cell  face  and  u^  is  an  inter- 
polated value  of  the  velocity  component  normal  to  the  cell  face 
and  represents  the  velocity  at  the  interface  at  the  end  of  the  time 
step.  For  example,  considering  transport  from  cell  k into 
the  cell  above  it,  ka,  the  cell  boundary  velocity,  v^  \ is 
defined  to  be  a simple  average  of  the  velocities  of  the  two 
cells.  However,  we  want  the  mass  transport  velocity,  v,  to  be 
the  cell  boundary  velocity  at  the  end  of  the  time  step.  At. 

This  transport  velocity  therefore  is  the  velocity  of  the  material 
located  a depth  Az  = vAt  into  the  donor  cell  k (see  Figure  2.3). 
However,  since  the  transport  velocity  v,  when  adjusted  by  the 
axial  velocity  gradient,  must  equal  the  cell  boundary  velocity, 
we  have 


V (1 


At)  = V 
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ka 


f a1 

Figure  2.3  --The  relationship  between  v'‘  , the  cell  boundary 

velocity  at  the  beginning  of  TPHASE,  and  v,  the 
transport  velocity. 


[i  > v.LKaj  - vjjO  At1 
Az  J . 

Similarly,  the  radial  mass  transport  velocity  between  the  donor 
cell  k and  its  neighbor  on  the  right,  kr,  is  defined  to  be 

j [u(k)  <•  u(kr)] 

“ ■ |~i  ^ ■'iua  At] 


Calculated  transport  masses  are  subtracted  from  the  donor 
cell  mass  and  added  to  the  acceptor  cell  mass,  and  the  cell 
momenta  and  internal  energy  are  updated  accordingly.  Since  all 
transport  terms  are  computed  using  the  post-HPHASE  quantities, 
a given  cell  is  updated  only  after  the  transport  terms  have 
been  calculated  for  all  four  of  its  boundaries. 


2. 1.1.1  Equation  of  Motion,  (2.11) 

The  TPHASE  contribution  from  the  equation  of  motion  (2.11) 
is 


Aj(mUj) 


The  axial  and  radial  momentum 
separately. 


(pu.Uj)n.dS  . 

increments  will  be  discussed 


I 


Effect  on  Axial  Motion 


At  each  of  the  four  cell  faces  comprising  the  surface 
of  the  cell,  the  specific  axial  momentun,  u^ , being  transported 
Is  the  axial  velocity  of  the  cell  from  which  the  mass  moves,  i.e.. 
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I 
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Uj  “ v(kcl) 

where  kd  is  the  index  of  the  donor  cell.  Since  this  velocity 
is  constant  over  the  cell  face  it  can  be  removed  from  the 
integral  sign  and  we  obtain 

A.j.(mv)  = v(kd) 

The  quantity  in  brackets  is  seen  to  be  the  mass  transported 
through  the  cell  face  as  calculated  in  Section  2. 2. 3.1.  Thus 
for  each  cell  face. 


^-At  J pu^n^dsj 


A.j.(mv)  = v[kd)  6m 

The  total  axial  momentum  increment  for  the  cell  is  the  sum 
of  the  increments  obtained  for  each  of  the  four  cell  faces. 


Effect  on  Radial  Motion 

The  specific  radial  momentum,  Uj , being  transported  is 
the  radial  velocity  of  the  cell  from  which  the  mass  moves, 

i.e.  , 


Uj  = u(kd) 

where  kd  is  the  index  of  the  donor  cell.  By  an  analysis 
similar  to  that  for  the  axial  momentum  increment,  the  radial 
momentum  increment  is 

A.p(mu)  = u(kd)  6m 

for  each  of  the  cell  faces.  The  total  radial  momentum  incre- 
ment for  the  cell  is  the  sum  of  the  increments  obtained  for 
each  of  the  four  cell  faces. 


i«rT  Tii-v;  Vk.  • '• 
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2 . 2 . 3 . 3 Ene 


luation 


The  expression  for  the  transport  of  energy  from  Eq.  (2.12) 


i^)  = J (pu^E^)  n.dS  . 


A,p(mE 


To  evaluate  this  integral,  the  transported  specific  energy^ E^^ 
across  each  cell  face  is  taken  to  be  that  of  the  donor  cell, 
kd,  i . e . , 


- EjCkd)  * \ [(u(kd))^  * (v(kd))^] 


and  the  total  energy  which  is  transported  across  a given  inter- 
face is  therefore  the  product  of  this  specific  energy  and  the 
associated  transport  mass,  which  was  computed  above  in  the 
continuity  equation,  i.e.  , 

d 

Aj(mEj)  » E^6m 

The  total  energy  increment  for  the  cell  is  the  sum  of  the 
increments  obtained  for  each  of  the  four  cell  faces.  The 
internal  energy  is  obtained  by  subtracting  the  new  "kinetic 
energy  (obtained  from  the  updated  mass  and  momenta)  from  the 
updated  total  energy.  The  partitioning  of  the  updated 
internal  energy  among  materials  in  multimaterial  cells  is 
discussed  in  Section  4.6.3. 


2.3  MATERIAL  MODEL 

As  is  well  known,  tho  conservation  equations  (Eqs.  2.1 
through  2.3),  which  consist  (in  two  dimensions)  of  four 
equations,  are  insufficient  to  allow  determination  of  the 
eight  unknowns  (p,  u,  v,  P,  E™,  s , s , s ) contained  in 
them  and  must  be  supplemented  by  equations  defining  the 
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material  model  to  be  used.  The  material  model  employed  in 
the  HELP  code  can  represent  inert  materials  or  high  explosives 
(HE)  and  provides  a fairly  sophisticated  HE  detonation 
routine  allowing  for  the  detonation  of  multiple  HE  packages  and 
for  detonation  around  inert  obstacles.  The  code  also  includes 
strength  effects,  modeling  them  in  the  elastic-plastic  regime, 
and  a simple  tensile  failure  model.  By  bypassing  the  strength 
effects  subroutine,  SPHASE,  the  code  solves  problems  in  com- 
pressible fluid  mechanics.  A more  complete  discussion  of  these 
points  follows. 


2.3.1  Equation  of  State 

The  HELP  code  employs  three  different  equations  of 
state:  the  Jones-Wilkins-Lee  (JWL)  equation  of  state 

high  explosives,  the  Tillotson  equation  of  state for 


for 

for  inert 


materials,  and  the  ideal  gas  equation  of  state.  The  array 
dimensions  of  the  HELP  code  allow  up  to  30  materials  in  any 
problem.  The  following  convention  is  used  to  relate  the  30 
material  numbers  to  one  of  the  three  equations  of  state: 
materials  1-19  are  assumed  to  use  the  Tillotson  equation  of 
state;  material  20  is  assumed  to  be  an  ideal  gas;  and  materials 
21-30  are  assumed  to  be  high  explosives  and  use  the  JWL  equation 
of  state. 


2. 3.1.1  High  Explosives 

Once  an  explosive  cell  is  detonated,  its  pressure  is 
computed  using  the  Jones-Wilkins-Lee  (JWL)  equation  of  state. ^ 
The  equation  is 


P(Ej,p)  = A 


(l  - 


+ B 1 - 


+ apEj 


L 
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where 


specific  internal  energy  in  ergs/g 
density  of  detonation  products  in  g/cm^ 
density  of  undetonated  explosive  in  g/cm^. 


The  constants  p^,  A,  B,  a,  a and  3 for  four  explosives  are 
listed  in  Table  2.1.  The  *onstants  fur  these  four  explosives 
are  defined  via  DATA  statements  in  EQST,  and  the  material 
numbers  of  these  explosives  are:  21-Comp  B;  22-TNT;  23-OCTOL; 

24-PBX  9404. 

The  JWL  equation  of  state  has  a pressure  maximum  which 
is  potentially  hazardous  to  the  pressure  iteration  (see 
Section  4.6.1).  The  pressure  iteration  adjusts  the  densities 
of  the  materials  in  a multimaterial  cell  to  achieve  pressure 
equilibrium.  During  this  iteration  process,  it  is  possible 
for  the  code  to  reference  the  JWL  equation  of  state  with  a 
density  corresponding  to  a pressure  beyond  this  maximum, 
thereby  preventing  convergence  of  the  iteration.  In  order  to 
preclude  this  possibility,  a variable,  RHOMAX,  is  defined 
via  DATA  statements  in  subroutine  CDT,  This  is  the  maximum 
value  of  density  for  each  HE  to  be  used  in  the  pressure 
iteration.  The  values  of  RHOMAX  for  the  four  specified 
explosives  are  listed  in  Table  2.1. 

The  sound  speed,  C,  in  the  detonated  explosive  is  de- 
fined by  the  equation 


f m 

- Pp' 


c = Vr  P/p  • 

The  value  of  r (defined  by  DATA  statements  in  subroutine  CDT) 
as  well  as  the  detonation  velocity,  D,  (defined  by  DATA  state- 
ments in  subroutine  DETIME)  and  detonation  energy,  E^ , (de- 
fined by  DATA  statements  in  ADDENG)  for  the  four  named  explo- 
sives are  listed  in  Table  2.1. 
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2. 3. 1.2  Inert  Materials 

The  equation  of  state  used  in  HELP,  for  initially  con- 
densed materials,  is  that  due  to  J.  H.  Tillotson, modified 
to  give  a smooth  transition  between  condensed  and  expanded 
states.  For  the  condensed  states,  i.e.,  when  p/p^  > 1,  or 
for  any  cold  states,  E^<Eg,  the  equation  has  the  form 


P = = 


a + 


V 


^ » 1 


Ejp  + Ay  + By' 


J 


where 


, ! 


Ej * specific  internal  energy  in  ergs/g 
p * material  density  in  g/cm^ 
n * p/pQ 

y " p/Pq  - 1 

For  expanded  hot  states,  i.e,,  when  p/p^  < 1 and  Ej > E' 
the  equation  of  state  has  the  form 


ff 

K'.  » 


aEjP  + 


l)E,P 


^ ^ 


+ Aye 


-B(Pq/p-1) 


LEoH 


-a(p^/p-l) 


A smooth  transition  between  the  condensed  and  expanded  states 

Is  insured  by  a transition  equation  for  the  intermediate  region 

defined  by  E < E < E"  and  p/p  < 1.  This  blended  portion  of 
sis  o 

the  equation  of  state  has  the  form 

„ CK,-E,)  Pj  ♦ (E'-Ej)  P, 
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In  these  equations,  a,  b,  a,  3,  E^,  E^,  E^,  A,  B,  and  are 
constants  for  the  particular  material.  The  values  of  these 
constants  for  nineteen  materials  are  given  in  Table  2.2. 

For  multimaterial  cells,  an  iteration  is  used  to  deter- 
mine the  cell  pressure  and  the  densities  of  the  various 
materials  within  the  cell.  This  iteration  is  discussed  in 
Section  4.6.1. 

2. 3. 1.3  Ideal  Gas 

As  noted  earlier,  material  number  20  is  assumed  to  be 
an  ideal  gas  and  the  pressure  is  calculated  from  the  ideal 
gas  equation  of  state 

P = (Y-1)  P Ej 

and  the  sound  speed  given  by 

C » 

with  p being  the  density,  Ej the  specific  internal  energy,  and 
Y the  ideal  gas  constant.  The  normal  density,  p^,  of  the 
ideal  gas  is  defined  in  a DATA  statement  to  be  that  of  air  at 
standard  temperature  and  pressure  (0.00129  g/cm^)  and  is  used 
only  for  editing  purposes.  Since  y is  stored  in  a non-dimensioned 
variable,  GAMMA,  the  code  has  the  capability  of  handling  only  one 
ideal  gas.  Note,  however,  that  the  JWL  equation  of  state 
reduces  to  the  ideal  gas  equation  of  state  if  A * B ■ 0 and 
a = Y-1*  Also,  in  order  for  a material  using  the  JWL  equation 
of  state  to  have  the  ideal  gas  sound  speed,  it  is  necessary 
to  have  T = y* 
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2.3.2  Clastic-Plastic  Constitutive  Relation 


Section  2.2.1  discussed  the  modification  of  the  material 
velocity  and  internal  energy  fields  due  to  the  effects  of  the 
stress  deviators  in  subroutine  SPHASE.  This  section  will  dis- 
cuss the  elastic-plastic  material  model  incorporated  in  the 
HELP  code  to  calculate  the  strass  deviators. 

Subroutine  SPHASE  has  two  major  tasks:  first,  it  computes 

instantaneous  strain  rate  deviators  and  updates  the  stress 
deviators,  and  second,  it  updates  the  cell  momenta  and  energy  to 
account  for  the  effects  of  the  stress  deviators  (as  discussed  in 
Section  2.2.1).  This  section  will  discuss  the  first  task  of 
SPHASE:  the  computation  of  the  strain  rate  deviators  and  the 

updating  of  the  stress  deviators,  including  a description  of 
the  von  Mises  yield  condition. 

2. 3.2.1  Strain  Rate  Deviators 

As  a first  step  in  computing  stress  deviators,  the 
instantaneous  strain  rate  deviators  are  computed  from  the 
velocity  field,  i.e.,  for  cylindrical  coordinates 


• 

e 

* V 

1 

(u  + V 

+ -) 

zz 

y 

O 

X y 

Sr  = S 

1 

■ T 

(u  + V 
X y 

^ -) 
X-' 

Sz  ■ (“y  * . 

where  the  x and  y subscripts  denote  partial  differentiation 
with  respect  to  those  variables.  Here  the  necessary  space 
derivatives  of  the  velocity  field  are  computed  (centered  at 
cell  centers)  in  a straightforward  manner,  i.e., 
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The  above  calculation  gives  only  the  strain  rate  deviator 
for  that  material  which  is  at  cell  center  at  the  beginning  of 
the  time  step,  and  it  is  necessary  to  recognize  that  constitutive 
equations  are  to  be  applied  only  along  a particle  path.  To 
this  end,  one  must  take  a convective  derivative.  Consider  the 
time  n position  q”  of  a particle  which,  at  time  n+1,  will  be 
at  a cell  center.  This  point  will  be  offset  from  the  cell 
center  by  6x  and  6y  where 


6x  ■ -u(k)  At 
6y  ■ -v(k)  At 


Ay 
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Now  we  wish  to  determine,  by  interpolation  using  the  value  of 
the  strain  rate  deviators  at  cell  centers,  the  strain  rate 
deviators  at  the  offset  point  q’^.  This  is  done  by  weighting 
the  contribution  of  neighboring  cells  in  proportion  to  their 
overlap  areas  with  a rectangle  of  cell  dimensions,  centered 
at  Q^.  For  example,  referring  to  the  above  sketch,  the 
weighting  factors  are 
for  cell  k 

I 

WK  *=  (Ax  - |6x|)  (Ay  - l6y|) 
for  cell  in  horizontal  direction 
WKH  = |6x|  (Ay  - |6y|) 


for  cell  in  vertical  direction 


WKV  = |6y  I (Ax  - |6x  j) 


for  cell  in  diagonal  direction 


WKD  = |6x|  |6y| 

so  that  the  resulting  interpolated  value  of  a strain  rate 
deviator  component,  say  is  given  by 


(WK)  ^^^(k)  + (WKH)e^^(kh)  + (WKV)E^^(kv)  + (WKD)ey^(kd) 


Having  thus  determined  , e'  , i'  for  the  mass  which  will 

zz  rr  rz 

be  at  the  center  of  cell  k at  the  end  of  the  time  step,  we 
can  now  calculate  the  strain  deviator  increments.  For 
example , 

Ae^_  = i'  At 
rz  rz 


2-31 


2. 3. 2.2  Stress  Deviators 


These  strain  deviator  increments,  ^re  next  used  in  the 

material  constitutive  equation  to  update  the  stress  deviator  tensor. 
To  do  this  one  needs  the  time  n stress  deviators  for  the 


=ij 


particle  in  question  and  these  s'^^j  are  determined  (from  stored 
cell  centered  values  of  the  s^^j  at  time  n)  by  an  identical 
area  weighting  procedure  as  described  above  for  strain  rate 


deviators, 


For  an  elastic-plastic  material  the  following  calculation 


then  (provisionally)  updates  the  s^j  according  to 


n+1 


ij 


■ =ij  * 2G  lej. 


where  G is  the  material  rigidity  modulus;  i.e.,  the  three 
components  of  interest  are  given  by 


s^'*'^  = s'  + 2G  Ae 
zz  z_  zz 


= s'  + 2G  Ac 
rr  rr  rr 


= s'^^  + 2G  Ae^^ 
rz  rz  rz 


The  provision  (above)  is  that  the  yield  condition  for  the 
material  not  be  violated. 
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2. 3. 2. 3 Yield  Criterion 


HELP  employs  the  von  Mises  yield  condition  which  requires 


that 


s.  .s.  . < 2Y‘- 
ij  11  - 

where  s s- • is  twice  the  second  invariant  of  the  stress  tensor^ 
ij  ij 

1 • C • f 

®lj  =ij  " =rr  ^ * *63  " 


where  s-, 
oO 


(s  + s ) is  the  hoop  stress  and  Y is  the 
'•  rr  zz 


shear  yield  strength.  If  the  yield  condition  is  violated 
during  a time  step,  the  deviatoric  stress  components  are 
proportionately  reduced  normal  to  the  yield  surface.  If  the 
yield  condition  is  not  violated,  the  deformation  is  purely 
elastic  and  the  stresses  are  not  reduced. 

A variable  yield  strength 


Y = t^o  ^1  ^ * ^2  ^ 5 - Ej/y 


is  defined  in  the  code  (in  subroutine  STRNG)  to  account  for 
the  increase  in  strength  at  high  compressions  and  the  decrease 
in  strength  at  elevated  energies.  In  this  equation  Y^, 

Yj^  and  Y2  are  input  constants,  p = p/p^  - 1 where  p is  the 
density  and  the  reference  density,  is  the  specific 
internal  energy  and  is  the  melt  energy.  This  expression 
is  never  allowed  to  be  negative,  however,  so  that  once  the 
melting  point  energy  is  exceeded  the  material  yield  strength 
is  set  to  zero.  It  siiould  oe  noted  tiiat  cells  (including  multi- 
material cells)  containing  material  which  has  failed  in 
tension  (see  the  following  subsection)  are  assumed  to  have  no 
shear  yield  strength.  Further  modifications  to  the  yield 
criterion  for  multimaterial  cells  are  discussed  in  Section  4.6.2, 
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2.3.3  Material  Failure  in  Tension 


The  tensile  failure  criterion  for  pure  cells  is  a simple 
one  based  on  relative  volume.  When  the  relative  volume 
(Vq/V  = p/Pq)  of  material  in  a pure  cell  reaches  a value 
which  is  greater  than  a specified  maximum  dilitation  (an 
input  con'stant  for  each  material,  AMDM-),  the  material  in 
the  cell  is  considered  to  have  failed.  A cell  which  has 
failed  is  not  allowed  to  sustain  tensile  hydrostatic  pressures 
nor,  as  noted  in  the  preceding  subsection,  to  have  shear  yield 
strength. 

The  tensile  failure  criterion  is  more  stringent  for  multi- 
material cells  in  that  they  are  not  allowed  to  sustain  negative 
pressures,  regardless  of  the  relative  volumes  of  the  materials 
they  contain,  since,  presumably,  the  negative  pressure  would  cause 
the  materials  to  separate  at  the  interface,  thus  immediately 
relieving  the  tension. 

2.3.4  Detonation  of  High  Explosives 


To  model  high  explosives  a detonation  procedure  has  been 
incorporated  into  HELP  as  well  as  special  prescriptions  which 

maintain  a sharp  boundary  between  the  detonated  and  undetonated 
explosive . 

2. 3. 4.1  Detonation  Procedure 

The  detonation  procedure  incorporated  into  HELP  makes  no 
assumptions  about  the  geometric  shape  of  the  explosive 
package,  the  number  of  explosives  in  the  problem,  or  the 
number  of  initiation  points,  although  certain  variable 
dimensions  allow  a maximum  of  ten  primary  and  ten  secondary 
initiation  points. 


2-34 


?v 

t 


Primary  initiation  points  may  be  located  in  any  con- 
figuration within  an  explosive  package,  and  may  have  a time 
delay  associated  with  them.  A primary  initiation  point  can 
be  detonated  at  any  time  on  or  after  time  zero  (the  start  of 
the  calculation).  For  example,  a string  of  five  detonators 
may  be  initiated  sequentially  or  simultaneously,  depending 
upon  the  desired  results. 

Secondary  initiation  points  which  are  used  to  detonate 
around  obstacles  may  also  be  located  in  any  configuration 
throughout  the  high  explosive  packages.  The  location  of 
secondary  initiation  points  should  be  selected  carefully  if 
they  are  to  perform  the  intended  function,  that  of  detonating 
around  obstacles.  The  initiation  times  for  each  of  the 
secondary  initiation  points  are  computed  by  the  same  method 
used  to  compute  the  detonation  time  for  a cell. 

Briefly,  the  numerical  method  for  handling  the  detonation 
of  a high  explosive  is  to  calculate  the  detonation  time  of 
each  explosive  cell  in  the  grid  based  upon  the  cell's  location 
and  the  constant  detonation  velocities  of  each  high  explosive 
through  which  the  detonation  front  passes  between  the  primary 
sources  and  the  cell  to  be  detonated. 

The  high  explosive  detonation  routine  employed  in  the 
HELP  code  allows  multi - initiation  of  multiple,  adjacent,  high 
explosive  packages.  For  the  case  of  adjacent  HE  packages 
the  current  method  is  restricted  to  explosives  whose  detonation 
velocities  and  C-J  pressures  are  within  ten  percent  of  one 
another  so  that  the  degree  to  which  one  of  the  detonations  is 
overdriven  is  minimal  and  can  be  neglected. 

Figure  2.4  illustrates  the  method  employed  to  compute 
the  detonation  times  of  cells  in  adjacent  high  explosive 
packages.  In  this  example  the  primary  initiation  point,  P^, 
is  in  the  HE.  package,  and  cell  K1  is  a pure  cell  in  the  HE, 


A,  <4^ 


package,  cell  K2  is  a pure  cell  in  the  IIE2  package,  and  cell 
K3  is  an  interface  cell  containing  both  HEj^  and  HE2.  In  all 
cases  the  time  to  detonate  a cell  is  the  sum  of  the  time  in- 
crements, t^,  associated  with  the  time  for  the  detonation 
front  to  travel  from  the  points  to  which  lie  on  the 

line  from  the  primary  initiation  point  to  the  center  of  the 
cell.  If  this  line  is  at  an  angle  which  is  less  than  or 
equal  to  45°  with  respect  to  tht'  x-rixis,  the  P.  arc  the  inter- 
cepts of  that  line  with  the  centers  of  the  grid  columns  (e.g., 
the  line  to  cell  K3) . If  the  angle  is  greater  tlian  45°,  the 
are  the  intercepts  of  that  line  with  the  centers  of  the 
grid  rows  (e.g.,  the  lines  to  cells  K1  and  K2)  . If  both  points, 
P^  j^  and  P^,  are  in  pure  HEj^  cells,  the  detonation  time,  t^, 
is  based  on  the  IlE^^  detonation  velocity,  Dj.  If  the  second 
point  P^,  is  in  an  interface  cell  containing  both  HEj^  and 
HE2  explosives,  the  detonation  time,  t^,  is  based  on  the 
maximum  of  the  HE^  and  ME2  detonation  velocities: 

^°1»  °2^max* 

Referring  to  Figure  2.4,  all  of  the  time  increments  for 
cell  K1  are  computed  using  . For  cell  K2,  however,  the  first 
three  increments  are  based  on  , the  next  two  on  (0i>^2^max’ 
and  the  last  one  on  D2.  For  cell  K3,  all  but  the  last  are 
based  on  and  the  last  is  a function  of  C^i>^2^max’ 

The  HELP  code  has  been  applied  to  the  jet  formation 
process  associated  with  lined  shaped  charges  in  which  the 
detonation  wave  has  been  "shaped"  by  the  insertion  of  inert 
material  into  the  high  explosive  package In  order  to 
accomplish  this  it  was  necessary  to  incorporate  the 
capability  of  detonating  around  obstacles  into  the  code. 

Figure  2.5  shows  a disc  of  inert  material  situated  in  a 
package  of  HE  in  which  detonation  is  initiated  at  a point  on 
the  axis  of  symmetry.  The  primary  initiation  point  is  shown 
as  Point  1.  The  detonation  time  of  each  cell  contained  in 
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area  I is  based  on  its  distance  from  Point  1 and  the  detona- 
tion velocity  of  the  HE.  Similarly,  the  cells  in  areas  II 
and  III  will  be  detonated  at  times  based  on  their  distances 
from  the  secondary  initiation  points,  2 and  3,  and  the 
detonation  velocity.  Points  2 and  3 in  Figure  2.5  act  as 
secondary  initiation  points  for  their  respective  areas  as 
soon  as  the  HE  in  the  cells  containing  these  points  has  been 
detonated . 

Both  the  primary  and  secondary  detonation  points,  as 
well  as  their  corresponding  detonation  regions,  are  part  of 
th.e  initial  conditions  and  are  specified  by  the  input  cards. 


At  the  beginning  of  the  problem  the  time  to  detonate  an 
explosive  cell  is  computed  and  stored  in  the  DETIM  array  (see 
Section  7.2.7).  On  the  cycle  for  which  the  value  of  T in  the 
calculation  equals  or  exceeds  the  detonation  time  associated 
with  a given  cell,  the  chemical  energy  of  detonation  of  the 
explosive  (E^  in  Table  2.1)  is  added  to  the  cell.  If  the  cell 
contains  more  than  one  explosive,  the  appropriate  chemical  energy 
of  detonation  of  each  high  explosive  is  added  to  the  cell. 

Each  cycle,  ADDENG  prints  the  row  (I)  and  column  (J) 
of  every  cell  detonated  on  that  cycle.  The  total  energy 
released  on  that  cycle  is  also  printed,  and  the  theoretical 
energy  of  the  grid  (ETK)  is  updated  accordingly. 


2. 3.4.2  Definition  of  the  Detonation  Front 

Because  the  detonation  procedure  used  by  HELP  detonates 
all  of  the  explosive  in  a cell  at  one  time,  the  detonation 
front  is  approximated  by  the  boundaries  between  detonated 
and  undetonated  cells.  In  order  to  prevent  small  numeri- 
cal signals  from  propagating  through  the  grid  ahead 
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of  the  detonation  front,  it  is  necessary  to  place  certain 
restrictions  on  cells  containing  undetonated  explosive  dur- 
ing the  three  phases  of  the  code. 

There  are  two  restrictions  necessary  to  prevent  the 
pressure  effects  from  propagating  undesired  numerical  sig- 
nals. The  first  of  these  is  implemented  in  CDT  and  requires 
that  any  cell  containing  undetonated  explosive  have  zero 
pressure.  For  these  cells,  the  equation  of  state  routine, 
EQST,  is  not  called  by  CDT.  This  restriction  is  necessary 
since  the  explosive  equation  of  state  applies  only  to  the 
detonated  explosive.  The  second  restriction,  implemented 
in  the  pressure  effects  subroutine,  HPHASE,  requires  the 
pressure  and  velocities  to  be  zero  at  the  four  boundaries 
of  an  interface  or  pure  cell  containing  undetonated  explosive, 
even  if  one  or  more  of  the  neighbor  cells  contains  detonated 
explosive. 

In  the  strength  effects  phase,  SPHASE,  any  mixed  or  pure 
cell  which  contains  explosive  is  assumed  to  have  no  strength 
and  its  deviator  stresses  are  set  to  zero.  This  assumption 
applies  to  cells  containing  detonated  as  well  as  undetonated 
explosive. 

In  subroutine  TPHASE,  the  transport  phase  of  the  code, 
no  mass,  momenta,  or  energy  is  transported  into  or  out  of  a 

cell  containing  undetonated  explosive.  These  transport 
restrictions  are  coded  into  DMCALC  as  well,  which 

computes  the  mass  of  each  material  to  be  transported  across 
the  four  boundaries  of  an  interface  cell.  If  an  interface 
cell  contains  undetonated  explosive,  no  mass  of  any  material 
is  transported  into  or  out  of  the  cell. 
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CHAPTER  III 


GRID  BOUNDARY  CONDITIONS 


3.1  BASIC  ASSUMPTIONS 

Additional  specifications  are  required  for  cells  which 
border  the  grid  boundaries  because  necessary  quantities  are 
not  defined  for  neighbor  cells  which  would  be  outside  the  grid. 
The  code  provides  for  a transmittive  or  reflective  boundary 
condition  at  the  bottom  and  a reflective  boundary  condition  at 
the  left  (see  Figure  3.1).  Boundary  conditions  for  grid  bound- 
ary cells  are  then  derived  by  assuming  fictitious  neighbor  cells 
outside  the  grid.  For  transmittive  boundaries  the  flow  variables 
are  the  same  in  the  fictional  cell  as  in  the  border  cell,  and 
for  reflective  boundaries  the  state  is  assumed  to  be  the  same  in 
the  fictional  cell  except  that  the  velocity  component  normal  to 
the  boundary  has  the  opposite  sign. 

These  assumptions,  as  they  apply  to  subroutines  ‘ oE, 
HPHASE,  and  TPHASE,  are  described  in  more  detail  in  the 
sections  that  follow.  Special  prescriptions  for  moving 
tracer  particles  on  and  near  the  grid  boundaries  as  well  as 
special  edited  quantities  which  indicate  the  work  done  and 
the  mass,  momentum,  and  energy  fluxes  at  the  grid  boundaries 
are  described  as  well. 


3.2  STRENGTH  PHASE  (SPHASE) 

The  SPHASE  grid  boundary  conditions  will  be  discussed 
as  they  pertain  to  the  three  subtasks  of  SPHASE: 

1.  the  calculation  of  cell  centered  strain  rates; 

2.  the  interpolation  of  cell  centered  stresses 
and  strain  rates  to  account  for  convection; 

3.  the  application  of  deviator  stresses  to  update 
cell  velocities  and  internal  energies. 
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Within  each  subtask  the  transmittive  and  reflective  boundary 
cases  will  be  discussed. 


3.2.1  Definition  of  Strain  Rate  Derivatives  for 
Cells  at  a Grid  boundary 


The  strain  rates  are  computed  using  space  derivatives 
centered  at  cell  centers.  For  example  the  finite  difference 
analogs  of  8v/9y  and  3u/3y  for  cell  k are 

3v  _ v(ka)  - v(kb) 

3y  “ -SAy.  ^ + Ay . + .SAy.^^ 


3u  ^ u(ka)  - u(kb1 

3y  “ .SAy..^  ♦ Ay^  + -SAy^^^  ’ 

with  ka  and  kb  being  the  indices  of  the  cells  above  and  below 
cell  k,  respectively,  and  Ay  being  the  axial  dimension  of  the 
corresponding  cell.  However,  if  cell  k is  at  the  top  trans- 
mittive grid  boundary,  these  space  derivatives  are  approxi- 
mated by 


3y  -5Ay^._^  + 


v(kl  - v(kb) 
Ay/,  + l.SAy, 


u(k)  - u(kb] 


•5Ayj_;^  + 


•SAy^ 


which  assumes  that  the  fictitious  cell  above  cell  k has  the 
same  velocity  components  and  axial  dimension  as  cell  k.  If, 
on  the  other  hand,  cell  k is  at  the  bottom  grid  boundary  and 
that  boundary  is  reflective,  these  space  derivatives  become 
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3v  _ y(ka)  r (-v(k)) 
8y  ■ l.  SAyV  + . SAy/^j 


3u  _ u(ka)  - u(k) 

3y  " l.SAy^  + .SAyV^^j 


which  assumes  that  a fictitious  cell  below  the  grid  has  the  same 
velocity  components  and  axial  dimension  as  cell  k,  except  that 
the  velocity  component  normal  to  that  boundary  has  the  opposite  sign 

Similar  approximations  are  made  for  8u/3x  and  3v/3x  for 
cells  at  the  right  transmittive  grid  boundary  and  for  cells 
at  the  left  reflective  grid  boundary. 


3.2.2  Definition  of  Interpolated  Strain  Rates  and 
Stresses  for  Cells  at  a Grid  Boundary 

SPHASE  computes  the  deviator  stresses  of  a particle  that 
will  be  at  the  center  of  cell  k at  the  end  of  the  time  step. 
The  strain  rates  and  time  n stresses  of  such  a particle  are 
computed  by  taking  an  area  weighted  average  of  cell  centered 
strain  rates  and  stresses.  (See  Section  2. 3. 2.1.)  If  cell 
k is  at  a grid  boundary,  and  the  velocity  field  is  such  that 
the  material  is  flowing  from  outside  the  grid,  then  it  is 
assumed  that  the  fictitious  cells  outside  the  grid  have  the 
same  stresses  and  strain  rates  as  the  contiguous  cells  in- 
side the  grid. 


3.2.3  Definition  of  Velocities  and  Deviator 


Stresses  at  Grid  Boundaries 


The  cell  velocities  and  internal  energies  are  updated 
in  SPHASE  using  cell  boundary  stresses  and  velocities  which 
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are  averages  of  the  cell  centered  quantities.  For  example, 
the  stresses  and  velocities  at  the  right  boundary  of  cell  k» 
which  is  not  at  a grid  boundary,  are 


-r  . u(r)  * u(kr: 


-r  . vfk)  * v(tr: 


_ s^^(k)  * s^^(kr) 


S “ ■ ' ' ' ' • 

rz  2 

If,  however,  cell  k is  at  the  right  transmittive  grid  boundary 
then  the  stresses  and  velocities  of  the  fictitious  cell  outside 
the  grid  are  assumed  to  be  the  same  as  those  of  cell  k,  and  the 
boundary  stresses  and  velocities  are 

r.r  _ ,,  flr-J 


u"  = u(k) 
v^  * v(k) 


^rr  ^rr^^^ 


And  if  cell  k is  at  the  left  reflective  grid  boundary, 
the  following  equations  define  the  stresses  and  velocities 
at  the  left  boundary: 


vCk) 


= 'rrC') 
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i.e.,  the  shear  stress  and  the  velocity  component  normal  to  the 
boundary  are  set  to  zero,  and  the  normal  stress  and  tangential 
velocity  component  are  assumed  to  be  equal  to  the  cell  k values. 
Likewise,  the  stresses  and  velocities  at  the  bottom  grid 
boundary,  when  it  is  reflective,  are 

u^  = u(k) 

.-.b  = 0 


= 0 
rz 


s„00 


3.2.4  Correction  to  the  Theoretical  Ener 


one  at 


for  Work 


The  theoretical  energy  of  the  grid,  ETH,  is  updated 
only  for  the  work  done  at  the  transmittive  grid  boundaries. 

At  the  right  boundary 

<5E^  = [s^^(k).u(k)  + s^^(k).vCk)  )-Aj.-At 

where  s^^,  s^^,  u,  v,  A^  are  the  deviator  stress  normal  to  the 
right  boundary,  the  shear  stress,  the  velocity  component  normal 
to  the  right  boundary,  the  velocity  component  tangent  to  the 
right  boundary,  and  the  area  of  the  right  cell  face,  respectively 
Similarly,  at  the  top  transmittive  grid  boundary 


^t  “ [s^^(k)*v(k)  + s^^(k)-u(k)]  *A^-At 


and  at  the  bottom  grid  boundary,  when  it  is  transmittive 


'vH-ff  / iwjjr  !■*»  *•*»'  '"•‘.t " ■ 


“ ts^^(k).v(k)  + s^j,(k) -u(k)]  .A^j.At 


When  the  bottom  boundary  is  reflective,  = 0,  since  in 

that  case  the  normal  velocity,  v,  and  the  shear  stress, 

s^„,  are  zero, 
rz 

3.3  HYDRODYNAMIC  PHASE  (HPHASEl 

The  effects  of  pressure  are  computed  by  using  pressures 
and  velocities  defined  at  cell  boundaries.  If  a cell  boundary 
is  part  of  the  grid  boundary,  the  definition  of  the  pressure 
and  velocity  at  that  boundary  depends  on  whether  or  not  that 
boundary  is  reflective  or  transmittive . 


3.3.1  Definition  of  Velocities  and  Pressures  at 
transmittive  Grid  Boundaries 

If  the  grid  boundary  is  transmittive,  the  cell  boundary 
pressure  is  assumed  to  be  equal  to  the  cell  centered  pressure. 
Likewise,  the  velocity  component  normal  to  the  boundary  is 
assumed  to  be  equal  to  the  cell  centered  velocity  component. 
For  example,  if  cell  k is  at  the  top  transmittive  boundary 


P(k) 


V®  = v(k) 


and  if  cell  k is  at  the  right  transmittive  boundary 


= P(k) 


u 


u(k) 


Because  the  transmittive  boundary  condition  assumes  all 
material  outside  the  grid  is  at  the  same  state  as  the  material 


in  the  border  cells  of  the  grid,  the  user  should,  by 
rezoning,  prevent  a strong  shock  from  moving  across  a trans- 
mittive  boundary.  In  such  an  instance  a numerical  signal  can 
be  generated  by  the  boundary  condition,  and  the  stresses  in 
the  material  behind  the  shock  can  become  noisy  and  meaningless. 


3.3.2  Definition  of  Velocities  and  Pressures  at 
deflective  tjrid  boundaries 

If  the  grid  boundary  is  reflective,  the  cell  boundary 
pressure  is  assumed  to  be  equal  to  the  cell  centered  pressure. 
However,  the  velocity  component  normal  to  the  boundary  is 
assumed  to  be  zero. 

For  example,  if  cell  k is  at  the  left  reflective  grid 
boundary  (an  axis  of  symmetry  in  cylindrical  coordinates) 

= P(k) 


and  if  cell  k is  at  the  bottom  grid  boundary  which  is  acting 
as  a reflective  boundary  (CVIS  = 0), 

P^  = P(k) 


3.3.3  Correction  to  Theoretical  Energy  for  Work 
Done  at  ffrid  boundaries  in  HpHAsIE 

The  theoretical  energy  of  the  grid,  ETH,  is  updated  in 
HPHASE  for  the  work  done  at  the  transmittive  grid  boundaries. 
At  the  top  transmittive  boundary 
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6E 


^ = -P  (k>v(k>A^-At  , 


and  at  the  right  transmittive  boundary 


*=  -P  (k>u(k>A*  At  , 


where  A^  and  A^  are  the  areas  of  the  top  and  right  cell  faces, 
respectively.  If  the  bottom  boundary  is  transmittive 


= P(k>v(k>AgAt  , 

and  if  it  is  reflective,  dE^^  * 0,  since 


-b 


0. 


3.4  TRANSPORT  PHASE  (TPHASE) 

Material  is  transported  through  the  grid  by  computing 
mass,  momentum  and  energy  fluxes  at  each  of  the  four  boundaries 
of  every  cell.  If  one  of  these  boundaries  is  a grid  boundary, 
these  transport  terms  are  determined  by  the  following  con- 
ventions for  transmittive  and  reflective  boundaries,  and  the 
theoretical  energy  is  adjusted  accordingly. 


3.4.1  Transport  of  Mass,  Momentum  and  Energy  Across 
Transmittive  Grid  Boundaries 

The  velocity  of  the  material  at  and  beyond  a grid  boundary 
is  assumed  to  be  the  same  as  that  of  the  boundary  cell.  The 
transport  velocity  used  to  compute  the  mass  transported  across 
a grid  boundary  is  therefore  the  cell  centered  velocity.  How- 
ever, if  the  cell  centered  velocity  component  normal  to  the 
grid  boundary  is  such  that  the  material  is  moving  into  the 
grid,  the  mass  transport  term  at  that  boundary  is  set 
to  zero.  For  example,  if  cell  k is  at  the  top  transmittive 
grid  boundary  and  v(k)  > 0,  then 
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6nia  = v(k) . At.  A^.  p(k) 

On  the  other  hand, if  v(k)  $ 0 then  = 0.  Likewise,  if 

a 

cell  k is  at  the  right  transmittive  boundary  and  u(k)  > 0, 
then 

6m^  = u(k) .At. A^. p(k) 

If  u(k)  £ 0,  however,  = 0. 

For  the  case  where  6m  > 0,  the  momentum  and  energy  trans- 
ported across  grid  boundaries  are  debtrmined  by  the  donor 
cell  method,  i.e . , 


6 (mu)  = 6m>u(k) 

6(mE)  = 6m.E(k) 

The  code  does  not  have  provisions  for  feeding  material 
into  the  grid  across  any  of  the  three  transmittive  grid 
boundaries;  however,  TPHASE  can  be  easily  modified  by  the 
knowledgeable  user  to  do  so. 


3.4.2  Change  of  Momentum  jgt  Cells  at  Reflective 
Grid  Boundaries  in  TPHASE 

Even  though  no  mass  is  transported  across  a reflective 
boundary,  the  momentum  of  a cell  at  a reflective  boundary 
is  modified  in  TPHASE. 

A useful  way  to  think  about  reflective  boundaries,  such 
as  the  axis  or  grid  boundaries,  is  to  take  them  to  be  planes 
(or  a line)  of  symmetry  through  which  equal  and  opposite 
transports  occur  in  both  directions. 
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Assume,  for  example,  that  the  mass  transport  calculation 
at  the  axis  indicates  that  a mass  6m*  would  pass  through  the 
axis  (i.e.,  u < 0) . Then  an  equal  and  opposite  transport  is 
assumed  from  the  left.  The  net  effect  on  cell  mass  is  zero. 

However,  there  is  an  effect  on  the  cell  momentum.  At  the 
axis,  momentum  6mu  is  lost  and  a momentum  6m(-u)  is  gained. 

The  net  result  is  a momentum  change  of  -26mu.  Thus  the  momentum 
equation  has  to  have  this  correction  term;  specifically 

A^mu  ^6  (m^^Uj^)  - 26mu  , 
i 

where  the  subscript  i denotes  the  four  sides  of  the  cell. 
Similarly,  it  can  be  shown  that  the  correction  term  at  the 
bottom  grid  boundary,  when  it  is  reflective,  is  -26mv. 

In  the  case  of  multimaterial  cells,  the  energy  equation 
must  also  be  modified.  For  example,  consider  a cell  on  the 
axis  of  symmetry,  where  u is  its  pre-TPHASE  radial  velocity 
and  u*  is  its  post-TPHASE  radial  velocity.  The  correction  to 
the  cell's  thermalized  kinetic  energy  (see  Section  4.6.3)  is 
as  follows: 


- ^ (u-u')^  + (-u-u')^  * 26muu' 

where  the  first  term  corresponds  to  the  mass  that  "left"  with 
velocity  u,  and  the  second  term  corresponds  to  the  mass  that 
"entered"  with  velocity  -u.  Similarly,  the  correction  term 
for  a multimaterial  cell  at  the  bottom  grid  boundary,  when  it 
is  reflective,  is  26mvv'. 


* The  discussion  above  is  independent  of  how  one  calculates  dm. 
However,  in  HELP,  at  the  axis  6m  ■ pu(iTArAz) At  if  u < 0, 
and  6m  ■ 0 if  u > TT7  The  area  term,  irArAz,  is  the  area  of 
the  surface  whicK  cuts  through  the  center  of  the  cell. 
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3.4.3  Correction  to  Theoretical  Energy  for  Energy 
l^ransported  Across  Grid  Boundaries  in  TPHASE 

In  TPHASE  the  theoretical  energy  of  the  grid,  ETH,  is 
adjusted  to  account  for  the  total  energy  of  the  mass  that 
leaves  the  grid  across  the  transmittive  grid  boundaries.  For 
example,  if  cell  k is  at  the  top  transmittive  grid  boundary, 
and  6m^  > 0,  then 

6Et  = 6m^-(E(k)  + 1/2  (u(k)^  + v(k)^)) 

where  E(k)  is  the  specific  internal  energy  of  cell  k.  There 
is  no  mass  loss  and  therefore  no  energy  loss  at  a reflective 
boundary. 


3.5  TRACER  PARTICLE  MOTION  NEAR  GRID  BOUNDARIES 

As  elsewhere  in  HELP,  subroutine  MOVTCR  assumes  the 
velocity  of  material  outside  a transmittive  grid  boundary 
is  equal  to  that  of  the  material  in  the  cells  bordering  the 
grid  boundary.  Therefore,  if  the  overlap  area  (described  in 
Section  4.1.2)  used  to  compute  the  velocity  of  a given  particle 
extends  into  a fictitious  cell  outside  the  grid,  the  velocity 
of  that  fictitious  cell  is  assumed  to  be  equal  to  its  neigh- 
bor inside  the  grid.  For  example,  as  illustrated  in  the 
following  figure,  the  radial  velocity  of  the  particle  is 
defined  as  follows: 

Aj^u(k)  + A2u(k)  + AjU(ka)  + A^u(ka) 

u = — — — 

Ai  + A2  + A3  + A4 
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I£  the  grid  boundary  is  reflective,  however,  the  velocity 
component  normal  to  that  boundary  in  the  fictitious  cell  is 
assumed  to  have  the  opposite  sign.  If  the  particle  is  near 
the  axis  of  symmetry,  as  illustrated  in  the  following  figure, 
the  radial  velocity  of  the  particle  is  defined  as  follows: 


A^u(k)  A2  (-u(k))  * Ag  (-u(ka))  + A^u(ka) 
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An  analogous  procedure  is  followed  for  computing  the  axial 
component  of  velocity  for  a particle  near  the  bottom  grid 
boundary  when  it  is  reflective  (CVIS  = 0) . 

The  coordinates  of  the  tracer  particles  are  stored  in 
grid  units.  If  the  new  position  of  a particle  is  computed 
to  be  beyond  the  grid,  its  x-coordinate  and/or  y-coordinate 
is  set  equal  to  the  grid  boundary  line  it  has  exceeded,  i.e., 
0,  IMAX,  or  JMAX,  since  it  is  not  possible  to  track  a 
particle  outside  the  grid.  Furthermore,  once  a particle  is 
moved  onto  a grid  boundary,  it  is  not  moved  off  of  that 
boundary,  only  along  it. 

3.6  EDITING  OF  FLUXES  AT  GRID  BOUNDARIES 

In  SPHASE  and  HPHASE  the  total  work  done  at  each  of  the 
transmittive  grid  boundaries  is  accumulated,  as  well  as  the 
internal  energy  lost  by  each  material  package  due  to  work 
done  at  these  boundaries  (calculated  in  subroutine  EOUT) . 

The  total  work  is  printed  by  subroutine  EDIT  under  the 
heading  WORK  DONE  and  is  included  in  the  total  printed 
under  the  heading  IE  OUT. 

In  TPUASE,  the  mass  that  leaves  the  grid  and  its  asso- 
ciated total  enfjrgy  and  momenta  are  accumulated  and  are 
printed  by  subroutine  EDIT  under  the  headings  MASS  OUT,  ENERGY 
OUT,  MU  OUT,  and  MV  OUT,  respectively.  The  kinetic  and 
internal  energy  lost  by  each  material  package  due  to  trans- 
port across  grid  boundaries  is  also  accumulated  and  is  in- 
cluded in  the  totals  printed  by  EDIT  under  the  headings 
KE  OUT  and  IE  OUT.  . 
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CHAPTER  IV 

MATERIAL  INTERFACES  AND  INTERFACE  CELLS 


The  following  sections  describe  the  HELP  method  of  de- 
fining and  propagating  material  interfaces  through  an 
Eulerian  grid  and  discuss  the  special  treatment  given  to  multi- 
material cells.  The  use  of  tracer  particles  to  define  the  in- 
terfaces and  the  manner  in  which  these  particles  are  moved  are 
discussed  in  Section  4.1.  The  method  for  storing  information 
on  several  materials  in  an  interface  cell  and  for  converting 
a pure  cell  into  an  interface  cell  is  described  in  Section  4.2. 
The  procedures  for  calculating  fractional  cell  face  areas 
(used  in  the  transport  of  materials  across  interface  cell 
boundaries)  are  given  in  Section  4,3.  The  adjustment  of  mass 
transport  terms  to  exactly  evacuate  a material  when  an  inter- 
face leaves  a cell  and  to  prevent  overemptying  of  a material 
is  discussed  in  Section  4.4.  In  Section  4.5  the  character- 
istics of  the  three  types  of  interface  cells  - multimaterial, 
free  surface  and  slipline  - are  given.  Finally,  the  special 
treatment  necessary  to  define  various  cell  and  material  prop- 
erties for  multimaterial  cells  is  given  in  Section  4,6. 


4.1  DEFINITION  OF  MATERIAL  INTERFACES 

4.1.1  Use  of  Tracer  Particles 

A set  of  massless  tracer  particles  are  initially  posi- 
tioned along  the  boundary  of  each  material  package.  These 
tracers  are  numbered  such  that  the  package  is  on  the  left  as 
one  proceeds  between  any  two  consecutive  tracers.  The  material 
interface  therefore  is  defined  by  the  line  segments  connecting 
these  tracers.  The  initial  spacing  of  the  tracers  is  specified 
when  they  are  generated  by  TSETUP  and  can  be  maintained  in  a 
specified  region  as  the  calculation  progresses  by  periodic 
calls  to  ADDTCR.  (See  description  of  ADDTCR  in  Section  8.2.) 
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To  prevent  the  boundaries  of  adjacent  packages  from  cros-  [ 

sing,  the  tracers  along  the  common  boundary  coincide  exactly.  j 

Because  these  identical  points  are  moved  with  exactly  the 

same  velocity,  they  remain  superimposed  during  the  entire  " 

calculation.  | 

To  allow  for  spatially  disconnected  subpackages,  the  j 

code  employs  the  following  convention.  The  last  point  of  a 
material  package  or  subpackage  is  a dummy  point,  its  x- 
coordinate  is  -1000,  its  y-coordinate  is  0.  A given  material 
package  can  be  divided  into  any  number  of  spatially  discon- 
nected subpackages. 


4.1.2  Movement  of  Tracer  Particles 

Each  material  package  is  circumscribed  by  a series  of 
massless  tracer  particles,  which  are  propagated  each  time 
step  a distance'^ 


Ax^  = u^  ^^INFACE 


(4.1) 


where  u.  is  a local,  density  weighted,  average  velocity  vector 
for  the  continuum,  determined  by  an  area  overlap  method  which- 
gives  weight  to  velocities  in  the  surrounding  cells.  Specifi- 
cally, a rectangle  of  cell  dimensions  is  superimposed  on  the 
particle  to  be  moved  and  then 


u, 


2 w(L) 


(4.2) 


Since  the  tracer  coordinates  are  in  grid  line  units,  this 
distance  is  converted  from  centimeters  to  cell  units. 
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where  w(L)  is  the  cell-centered  velocity*  of  the  overlapped 
cell  L,  is  the  overlap  area  and  Pj^  is  the  total  cell  mass 
divided  by  the  total  cell  volume.  This  procedure  gives  a 
spatially  continuous  velocity  field  for  particle  propagation. 

The  density  weighting  causes  the  particles  to  follow 
the  motion  of  the  denser  material  and  prevents  the  inter- 
face from  becoming  unstable  in  regions  where  there  are  large 
density  discontinuities. 

Infrequently,  the  sum  of  the  area  terms,  A^,  equals  zero, 
which  indicates  all  of  the  overlap  cells  associated  with  a 
given  tracer  particle  are  empty  and  it  cannot  be  moved  in  the 
usual  manner.  When  this  occurs,  the  code  will  move  the 
isolated  tracer  to  the  position  of  one  of  its  immediate 
neighbors.  However,  if  the  isolated  tracer  is  on  a grid  bound- 
ary, it  will  not  be  moved  off  of  that  boundary,  only  along  it. 


If  INFACE  is  subcycled  (CYCMX  > 1),  then  the  time  step. 


At 


INFACE' 


used  to  move  the  tracers  on  a given  subcycle,  is  the 


time  step  computed  by  CDT  divided  by  the  number  of  INFACE 
subcycles,  i.e. 


At 


At 


CDT 


INFACE 


CYCMX 


(4.3) 


*These  cell-centered  velocities  have  been  updated  by  HPHASE 
(pressure  effects)  and  SPHASE  (strength  effects)  for  this 
time  step. 
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The  passive  tracers  (XP,  YP) , used  only  for  tracking 
the  mass  within  the  package  boundaries,  are  moved  by  the 
same  method,  except  that  if  INFACE  is  subcycled  they  are 
moved  only  on  the  last  subcycle  using  the  time  step 


4.2  CREATION  OF  INTERFACE  CELLS 

In  the  discussion  that  follows  an  interface  cell  is  any 
cell  whose  boundary  is  intersected  by  a material  interface, 
whereas  a pure  cell  is  completely  within  one  of  the  material 
package  boundaries.  An  interface  cell  may  contain  one  or 
more  distinct  materials,  it  may  contain  a void  and  therefore 
be  a free  surface  cell,  and  it  may  contain  a slipline.  The 
characteristics  of  the  various  types  of  interface  cells  are 
given  in  Section  4.5.  The  way  in  which  interface  cells  are 
flagged  and  the  additional  information  stored  for  each  of 
them  are  discussed  in  the  sections  below. 
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4.2,1  Use  of  MFLAG  Array 

In  all  cases,  the  value  of  MFLAG (k)  of  an  interface  cell 
k is  greater  than  100.  The  value  of  MFLAG  (k)  of  a pure  cell 
k,  however,  is  less  than  or  equal  to  the  number  of  material 
packages  in  the  grid  (0,  1,  ...,  NMAT) . 

For  interface  cells,  the  MFLAG  value  has  a twofold 
function.  First,  as  indicated  above,  it  indicates  that  the 
cell  is  intersected  by  an  interface.  Second,  it  links  the 
k index  of  the  cell  to  an  index  in  those  arrays  which  store 
the  mass,  density,  specific  internal  energy,  and  velocity 
components  of  each  material  in  the  interface  cell. 

Given  m = MFLAG(k)  - 100,  the  mass  of  material  package 
n in  cell  k is  stored  in  XMASS  (n,m) ; similarly  the  density, 
specific  internal  energy,  and  velocity  components  of  material 
package  n are  stored  in  RIIO  (n,m),  SIE  (n,m),  US  (n,m)  and 
VS  (n,m),  respectively.  In  this  way  the  k index  of  the  cell 
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arrays  are  linked  to  the  m index  of  the  material  arrays. 
Furthermore,  the  cell  arrays  are  defined  by  summing  over  the 
materials  in  the  interface  cell  as  follows: 


AMX(k)  XMASS(n,m) 


U(k)  = yuS(n,m)‘XMASS( 


/ n 


(4.4) 


XMASS(n,m)  (4.5) 


V(k)  = J]vS(n,m)-XMASS(n,m)/^  XMASS(n,m)  (4.6) 

n / n 

AlX(k)  = ^SIE(n,m)  •XMASS(n,m)/j^XMASS(n,m).  C4.7) 

n / n 


Listed  in  the  table  below  are  the  arrays  which  store  cell 
quantities  and  the  arrays  which  store  material  quantities 
within  the  interface  cells.  Throughout  this  document  these 
will  be  referred  to  as  cell  arrays  and  material  arrays,  re- 
spectively. 


QUANTITY 

CELL  ARRAY 

MATERIAL  ARRAY 

Mass 

AMX 

XMASS 

Specific  Internal  Energy 

AIX 

SIE 

Radial  Velocity 

U 

US 

Axial  Velocity 

V 

VS 

Density 

(Not  stored) 

RHO 

Pressure 

P 

(Not 

stored) 

Shear  Stress 

STRSRZ 

(Not 

stored) 

Radial  Stress  Deviator 

STRSRR 

(Not 

stored) 

Axial  Stress  Deviator 

STRSZZ 

(Not 

stored) 

4.2.2  Definition  of  Material  Arrays 

Three  subroutines,  NEWFLG,  NEWMIX,  and  NEWRHO  are  called 
to  define  the  material  arrays  for  a cell  k that  is  being  con- 
verted from  a pure  to  an  interface  cell,  which  occurs  when  a 
material  interface  initially  enters  a pure  cell. 

First  NEWFLG  is  called  to  find  an  unused  location,  m, 
in  the  material  arrays.  The  convention  which  denotes  that  a 
material  array  location  m is  not  being  used  is  that  RH0(l,m)=-l. 
Having  found  an  unused  location  (HELP  performs  an  error  exit 
if  no  such  location  is  found),  NEWFLG  redefines  the  value  of 
MFLAG(k)  to  be  m + 100. 

In  the  case  when  an  interface  cell  becomes  a pure 
cell,  which  occurs  when  all  material  interfaces  leave  the  cell, 
it  is  necessary  to  release  the  members  of  the  material  arrays 
which  have  stored  the  material  properties  for  that  interface 
cell.  This  process,  which  is  implemented  in  TPHASE,  consists 
of  initializing  the  particular  material  variables  to  zero  and 
setting  RH0(l,m)  = - 1,  which  denotes  that  the  material 
variables  whose  second  index  is  m are  unused  and  can  be 
assigned  to  another  interface  cell  by  NEWFLG.  The  process 
of  converting  an  interface  cell  to  a pure  cell  is  completed 
by  setting  the  cell's  MFLAG  value  to  the  package  number  of 
the  material  that  now  fills  the  cell,  or  to  zero  if  the  cell 
has  become  completely  void. 

NEWMIX,  after  calling  NEWFLG,  defines  the  material  arrays 
for  the  material  which  is  currently  in  the  cell.  If  mo  was 
the  value  of  MFLAG  (k)  on  the  previous  cycle, when  cell  k was 
pure,  and  m is  the  location  in  the  material  arrays  assigned 
to  cell  k by  NEWFLG  [m  + 100  = MFLAG(k)],  then  NEIVMIX  defines 
the  material  variables  for  the  material  in  the  cell  as 
follows  (providing  cell  k was  not  a pure  void  cell,  i.e., 
mo  f 0) : 
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XMASS  (mo,  m)  = AMX(k) 

SIE  (mo,  m)  = AlX(k) 

US  (mo,  m)  = U(k) 

VS  (mo,  m)  - V(k) 

RHO  (mo,  m)  = AMX(k)yCELL  VOLUME  . 

If  mo  = 0,  the  cell  was  previously  a void  cell  and  has  be- 
come a free  surface  cell,  since  it  contains  a material 
interface.  In  that  case  RHO  (NVOID,  m)  is  set  equal  to  1, 
signifying  that  the  cell  is  a free  surface  interface  cell. 

Any  interface,  of  course,  is  a boundary  between  two 
materials,  or  one  material  and  a void.  Therefore,  when  a 
cell  becomes  an  interface  cell,  both  materials  or  the  void 
have  to  be  identified.  The  identity  of  one  material,  or  void, 
is  indicated  by  the  flag  of  the  cell  when  it  was  pure  (mo  in 
the  example  above) . The  other  material  or  void  is  identified 
when  the  code  senses  which  package  tracer  particles  are  in- 
tersecting the  cell's  boundaries.  If  that  other  package  is 
the  void,  then  RH0(NV0ID,m)  is  set  to  1.  Otherwise  sub- 
routine NEWRHO  is  called  to  assign  a density  to  the  other 
material.  NEWRHO  picks  the  density  from  the  neighbor  cells 
nearest  the  intercept  of  the  interface  with  the  cell  boundary. 
This  density  value  is  used,  in  most  cases,  as  a flag  in- 
dicating which  other  material  interface  has  entered  the  cell. 
Since  the  densities  of  multimaterial  interface  cells  are  re- 
defined in  the  pressure  iteration  (see  Section  4.6.1)  at  the 
beginning  of  the  next  cycle,  the  value  chosen  need  only  be 
an  approximation  of  the  density  of  the  material  about  to  be 
transported  into  the  cell.  Only  if  the  interface  cell  had 
been  a pure  void  cell  (mo  = 0)  and  will  contain  only  the  new 
material  which  is  entering  the  cell  does  the  density  defined 
by  NEWRHO  remain  unmodified  by  the  pressure  iteration.  This 
is  the  only  case  in  which  the  density  assigned  to  the 
material  by  NEWRHO  is  used  as  a real  measure  of  the  density 
of  the  material  in  the  interface  cell. 
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It  is  important  to  note  that  the  material  density  array, 
RHO,  is  used  as  a flag  as  well  as  a material  property. 
RH0(n,m)>0  indicates  that  the  interface  of  material  package 
n intersects  the  m interface  cell.  This  flag  must  be  set 
before  any  mass  of  package  n can  be  transported  into  inter- 
face cell  m.  However,  it  is  possible  for  the  density 
variable,  RIIO  (n,m),  to  be  non-zero  and  the  corresponding 
mass  variable,  XMASS  (n,m),  to  be  zero.  In  that  case,  the 
non-zero  density  only  indicates  that  the  interface  of 
material  package  n is  intersecting  the  cell. 


On  the  other  hand,  when  the  material  n interface  leaves 
the  m interface  cell,  the  density  variable,  RHO  (n,m),  is  set 
to  zero  to  indicate  that  all  the  mass  of  material  n should 
be  "evacuated"  from  the  cell  on  that  cycle.  In  this  case 
the  density  variable,  RHO  (n,m),  being  zero  and  the  corres- 
ponding mass  variable,  XMASS(n,m) , being  nonzero  signifies 
that  any  influxes  of  that  material  into  the  cell  should  be 
set  to  zero  and  the  outfluxes  adjusted  so  that  their  sum  ex- 
actly equals  the  remaining  mass.  This  adjustment  of  the  trans- 
port terms  is  done  by  subroutine  DMADJ  and  is  discussed  in 
more  detail  in  Section  4.4. 


4.2.3  Accounting  for  Subcycles 

INFACE  usually  is  subcycled  to  minimize  the  transport  noise 
which  occurs  when  an  interface  leaves  a cell.  The  routine  is 
executed  two  or  more  times  each  cycle  using  a fraction  of  the 
time  step  on  each  subcycle.  It  is  therefore  possible  to 
create  an  interface  cell  after  one  or  more  subcycles  of 
INFACE  have  been  completed.  In  that  case  the  mass  transport 
variables  (SAMMP,  SAMMY,  SGAMC,  SAMPY)  for  that  cell  have 
not  been  accumulated  for  the  completed  subcycles  and  they 
need  to  be  computed  before  adding  in  the  transport 
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terms  for  the  current  and  subsequent  subcycles.  Therefore, 
when  INFACE  has  completed  at  least  one  subcycle  and  storage 
for  the  new  interface  cell  has  been  set  up,  NEWMIX  calls 
DMCALC,  where  the  computation  of  transport  terms  is  performed. 


4.3  CALCULATION  OF  CELL  FACE  AREAS  FOR  TRANSPI 
INTERFACE  CELL  BOUNDARIES 


The  computation  of  fractional  cell  boundary  areas  is 
performed  in  order  to  use  the  material  interface  positions 
to  define  the  transport  of  a given  material  across  an  inter- 
face cell  boundary.  For  pure  cells  the  mass  transport  be- 
tween two  cells  is 


V A At, 


(4.8) 


where  is  the  donor  cell  density,  v is  an  average  of  the 
cell-centered  velocities  of  the  two  cells,  A is  the  area  of 
the  cell  boundary  across  which  the  mass  is  being  transported 


and  Atpj,™  is  the  time  step*  computed  by  CDT. 


For  interface  cells  this  equation  becomes 


''"'n  ■ On  % % “iNFACE 


(4.9) 


where  the  density,  p^,  and  velocity,  v^^,  are  based  on  material 


quantities  (RHO,  US,  VS)  rather  than  cell  quantities  (AMX,  U, 


V) . The  area  term,  A , is  associated  with  material  n and  is 

n 


determined  by  the  position  of  the  material  n interface,  while  the 
time  step,  is  the  time  step  for  the  cycle  computed  by 

CDT,  divided  by  the  number  of  passes  through  subroutine  INFACE 
(see  Eq.4.3).  Thus  the  code  computes  a mass  transport  term  for 
each  material  which  intersects  any  boundary  of  an  interface  cell. 


* See  DT  in  Section  10.3, 
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This  section  describes  the  procedures  used  to  compute  the 
cell  face  area  terms  which  are  used  to  transport  material 
across  the  boundaries  of  interface  cells. 

4.3.1  Defining  Fractional  Areas  of  Intersected  Cell  Faces 

The  cell  face  area  terms  associated  with  each  material 
in  an  interface  cell  are  computed  in  subroutine  FRACS.  The 
method  given  below  for  computing  these  areas  puts  no  restrictions 
on  how  many  interfaces  are  in  a cell  nor  on  how  many  times  a 
single  interface  crosses  a cell  boundary.  It  does,  however, 
assume  that  a material  interface  never  intersects  itself, 
and  experience  has  shown  that  the  logic  of  the  code  breaks 
down  if  it  does.  This  is  not  a limitation  of  the  method 
since  the  boundary  will  not  cross  itself  unless  the  tracers 
become  too  sparse,  and  this  can  be  prevented  by  periodi- 
cally adding  new  tracers  in  subroutine  ADDTCR.  (See  Section  8.2) 

The  method  for  computing  the  area  terms  will  be  better 
understood  if  the  following  assumptions  and  conventions  are 
kept  in  mind.  First  of  all,  FRACS  computes  and  stores  only 
the  area  terms  for  the  right  and  top  boundaries  of  each  in- 
terface cell, since  the  cells  below  and  on  the  left  provide 
the  area  terms  for  the  other  two  boundaries.  FRACS  stores 
these  area  terms  for  interface  cells  in  two  arrays,  FRACRT  and 
FRACTP,  for  the  right  and  top  cell  boundaries,  respectively. 

These  arrays  are  doubly  dimensioned^  given  FRACRT (n,m)  or 
FRACTP(n,m),  the  first  dimension,  n,  gives  the  material  pack- 
age number,  and  the  second  dimension,  m,  is  the  index  that 
links  the  interface  cell  to  the  material  arrays  by  the 
relation  m = MFLAG(k)  - 100.  Also,  because  the  area 
term  of  a given  material  package  appears  in  the  mass  trans- 
port equation,  it  must  be  defined  for  the  right  and  top 
boundaries  of  every  interface  cell,  even  for  those  boundaries 
not  intersected  by  a material  interface. 


4-10 


Let  us  first  consider  how  FRACS  computes  the  area  terms 
for  the  cell  boundaries  that  are  intersected  by  a material 
interface.  The  set  of  tracer  particles  which  circumscribe 
a given  material  are  considered  consecutively  (two  at  a time) 
by  FRACS.  (TX  and  TY  are  the  FORTRAN  variable  names  for  the 
particles'  x and  y coordinates,  respectively.  The  TX  and  TY 
arrays  are  doubly  dimensioned,;  the  first  dimension  gives  the  mat- 
erial package  number,  the  second  gives  the  sequential  order- 
ing of  the  particles.  For  simplicity  the  coordinates  of 
the  tracers  are  in  grid  line  units  rather  than  centimeters, 
since  in  grid  line  units  the  integer  part  of  the  coordinates  in- 
dicates whether  a pair  of  tracers  straddles  one  or  more  cell 
boundaries.)  If  a pair  of  material  n tracers  straddles  a cell  bound 
ary,  FRACS  computes  the  intercept  of  that  boundary  with  the  line 
between  the  two  tracers.  Given  that  point  of  intersection 
FRACS  defines  the  area  term  for  material  n at  that  cell 
boundary  to  be  the  fractional  cell  boundary  area  on  one  side 
of  the  intercept.  To  determine  which  side  of  the  intercept 
to  associate  with  material  n, FRACS  uses  the  convention  that 
material  n is  on  the  left  as  one  proceeds  from  the  first  to 
the  second  tracer  of  the  pair. 

For  the  cell  boundary  not  intersected  by  a material 
interface  the  area  term  for  a given  material  is  either  the 
total  cell  boundary  area  or  is  zero.  The  following  discus- 
sion illustrates  how  the  area  terms  of  these  non- intersected 
cell  boundaries  are  correctly  defined  by  "presetting"  them 
when  an  interface  enters  a cell  and  "resetting"  them  when  it 
leaves. 


4.3.2  Presetting  Non- Intersected  Cell  Face  Areas 


The  sequential  order  of  the  tracer  particles  determines 
whether  an  interface  is  entering  a cell  or  leaving  it.  When  the 


interface  is  entering  a cell,  and  is  crossing  this  cell 
boundary  for  the  first  time,  the  program  processes  the  other 
boundaries  of  the  cell  in  a clockwise  order,  presetting  the 
corresponding  area  terms  to  the  total  area  of  that  cell  bound- 
ary, until  it  encounters  a side  that  has  already  been  crossed 
by  the  same  material  interface  (i.e.,  a side  which  has  an 
associated  area  term  that  is  non- zero  and  is  less  than  the 


total  area).  The  first  three  c 
the  presetting  procedure. 

Case  (1) 


ases  diagrammed  below  illustrate 

Area  terms  preset: 

Bottom  = IT  - X? 

. (y.  - y._j)  zn 

Top  » ifx?  - Xj_j] 

Area  terms  preset: 

FottOm  = TT  |x^  - 


i-1 


1 


Case  (3) 


Area  terms  preset: 

None  - the  bottom  boundary  has 
been  intersected  previously  by 
the  same  interface. 


In  order  that  there  be  no  restriction  on  the  number  of 
times  the  material  interface  crosses  a given  cell  boundary, 
the  presetting  procedure  is  generalized  in  the  following 
manner.  If  the  interface  of  material  n crosses  a cell 
boundary  for  the  second  time,  the  code  does  not  preset  the 
area  terms  of  the  other  sides  of  the  cell  if  material  n lies 
between  the  two  points  of  intersection.  This  is  illustrated 
by  Case  4. 

Case  4 


Area  terms  preset: 


None  - material  n 
lies  between  the 
two  intersection 
points. 


The  pro  , n senses  that  material  n is  between  the  points 
when  1 of  the  area  term  computed  on  the  previous 

cross,  as)  and  the  one  coinputed  currently  is  greater  than 
the  total  area  of  thn  cell  face. 

If,  as  in  Case  5,  the  Kuni  of  the  area  term  computed  on 
the  previous  crossing (s)  '’nd  the  one  computed  currently  is 
less  than  the  total  area,  material  n is  outside  the  two 
intercepts  and  the  procedure  described  above  for  presetting 
the  area  terms  of  the  other  sides  of  the  cell  is  followed. 


It  should  be  noted  that  these  conventions  apply  re- 
gardless of  how  many  times  the  interface  crosses  the  boundary, 
as  illustrated  by  Case  (6)  and  Case  (7). 


Case  (6) 


Area  terms  preset: 

/ 2 2 \ 

Bottom  = 

Left  = (yj  - 


Case  (7) 


f/XrC 


di: 


Previous 


- X - — Current 


Area  terms  preset: 


None 


vrr. 


t 


4.3.3  Resetting  Non- Intersected  Cell  Face  Areas 

The  procedure  for  resetting  the  area  terms  when  the  in- 
terface leaves  a cell  is  the  following.  If  the  interface  is 
crossing  the  cell  boundary  for  the  first  time,  the  program 
processes  the  other  faces  of  the  cell  in  a clockwise  order 
and  resets  their  area  terms  to  zero,  until  it  encounters  a 
side  that  has  been  crossed  previously  by  that  same  material 
interface.  This  procedure  is  illustrated  by  Cases  (8),  (9), 
and  (10)  below. 

Case  (8) 


Area  terms  reset; 


Top  ■ 0 


Right  = 0 
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Case  (9) 


Area  terms  reset; 


Top  = 0 


Casa  (10) 


Area  terms  reset 


None 


In  Cases  (11)  and  (12)  below,  the  side  where  the  inter- 
face leaves  has  been  previously  crossed  by  this  interface. 
The  program  sums  the  area  terms  resulting  from  any  previous 
crossings  of  that  boundary  and  the  current  crossing.  If,  as 
in  Case  (11),  the  sum  of  these  areas  is  less  than  the  total 
cell  boundary  area,  the  material  is  outside  the  two  inter- 
cepts and  none  of  the  area  terms  of  the  other  sides  are  set 
to  zero. 
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Case  (11) 


Area  terms  reset; 
None 


, Previous 


ft#r/al  n 


'Current 


If,  however,  the  sum  is  greater,  as  in  Case  (12),  the  material 
is  between  the  intercepts  and  the  program  proceeds  to  reset 
to  zero  the  areas  of  the  other  non- intersected  sides  of  the 
cell . 


Case  (12) 


Previous  Area  terms  reset; 


Current  xop  = 0 


Right  = 0 


Bottom  = 0 


4.4  ADJUSTMENT  OF  INTERFACE  CELL  MASS 


The  four  maS  transport  terms  for  interface  ceils  (one  for 
each  cell  face)  We  computed  by  DMCALC  using  material  densities 
and  velocities  and  the  cell  face  area  terms  associated  with  the 
material  packages  as  computed  by  subroutine  FRACS.  However,  sub- 
sequent to  this  calculation,  DMADJ  is  called  to  (possibly) 
modify  these  transport  terms.  There  are  two  reasons  the 
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transport  terms  may  need  to  be  adjusted:  one,  to  exactly  empty 

a cell  of  a material  when  its  interface  leaves  the  cell,  i.e., 
evacuation;  and  two,  to  prevent  the  outflux  of  material  from 
exceeding  the  mass  of  the  material  in  the  cell,  i.e.,  over- 
emptying. 


4.4.1  Evacuation  of  a Material 


To  signify  that  the  material  n interface  has  passed  out 
of  cell  k the  density  of  material  n,  RMOCn,m),  is  set  to  zero. 

This  zero  density  is  used  as  a signal  to  subroutine  DMADJ  to 
adjust  the  transport  terms  (one  for  each  cell  face)  of 
material  n so  that  cell  k will  be  exactly  evacuated  of 
material  n.  This  adjustment  of  the  four  transport  terms 
usually  occurs  twice  in  the  computational  cycle;  once  at 
the  beginning  of  the  first  INFACE  subcycle  (in  case  the 
material  r interface  left  the  cell  on  the  previous  cycle*)  and 
once  after  all  INFACE  subcycles  are  completed  (in  case  the  mate- 
rial n interface  moved  out  of  cell  k during  or  after  the  first, 
and  before  the  last,  INFACE  subcycle).  For  the  first  case  the  four 
transport  terms  from  the  previous  cycle  are  used  to  indicate  the 
direction  of  the  flow  and  the  direction  of  the  evacuation.  In  the 
second  case,  the  four  transport  terms  just  computed  by  DMCALC  are 
used  to  determine  how  material  n is  to  be  evacuated. 


If  material  n is  to  be  evacuated  from  cell  k,  then  all 
influx  transport  terms  of  material  n are  set  to  zero,  and 
each  outflux  term  is  adjusted  in  proportion  to  its  fraction 
of  the  total  outflux  computed  by  DMCALC.  If  the  total  out- 
flux computed  by  DMCALC  is  zero,  the  material  is  removed 
from  the  grid  (evaporated)  at  the  end  of  TPHASE. 


*The  tracer  particles  are  moved  at  the  end  of  each  INFACE  sub- 
cycle. If  the  interface  moves  out  of  a cell  on  the  last  sub- 
cycle of  INFACE,  the  code  does  not  sense  this  until  the  first 
pass  through  INFACE  on  the  next  computational  cycle. 
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These  evacuation  procedures  are  used  for  a cell  that 
has  lost  one  interface,  but  is  still  cut  by  another  inter- 
face, as  well  as  for  cells  that  have  become  pure.  After 
the  area  terms  associated  with  each  material  package  have 
been  computed,  FLGSET  searches  for  cells  that  were  interface 
cells  on  the  previous  cycle  but  are  no  longer  cut  by  an 
interface.  Such  a cell  has  become  pure  and  its  flag, 
MFLAG(k),  is  made  negative  to  signify  that  fact  until  the 
transport  has  been  completed  in  TPHASE.  To  determine  what 
material  package  a newly  pure  cell  k belongs  to, FLGSET 
first  looks  for  a neighbor  that  is  pure  and  assumes  cell  k 
will  belong  to  the  same  package  as  one  of  its  pure  neighbors. 
If,  however,  all  of  its  neighbors  are  interface  cells,  the 
top  cell  face  area  terms  of  the  cell  below?  kb,  are  examined. 
If  the  cell  face  area  for  material  n at  the  top  of  cell  kb 
is  equal  to  the  total  cell  face  area,  then  cell  k must  be  in- 
side the  material  package  n boundary.  Therefore,  all 
material  densities  for  cell  k except  that  for  material  n 
are  set  to  zero,  and  at  the  end  of  TPIIASE  cell  k becomes  a 
pure  material  package  n cell  with  MFLAG(k)  = n. 

4.4.2  Overemptying  of  a Material 

After  the  subcycles  of  INFACE  are  completed,  subroutine 
DMADJ  is  called  to  check  that  the  outflux  of  each  material 
in  each  interface  cell  does  not  exceed  its  current  mass. 

If  the  sum  of  the  outgoing  transport  terms  of  material  n is 
greater  than  the  current  mass  of  material  n,  the  outgoing 
transport  terms  of  material  n are  reduced  proportionately, 
as  in  the  evacuation  procedure  described  above. 


*The  cell  on  the  left  of  cell  k is  used  analogously  if  cell 
k is  in  the  bottom  row  of  the  grid. 
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4.5 


There  are  several  types  of  interface  cells.  In  the 
following  paragraphs  they  will  be  defined  and  then  discussed 
in  terms  of  special  prescriptions  that  apply  to  them. 

4.5.1  Multimaterial  Cells 

Most  interface  cells  contain  more  than  one  material,  and 
the  only  practical  restriction  on  the  number  of  materials  in 
any  one  cell  is  the  array  dimensions  imposed  by  the  user?  A 
detailed  description  of  the  determination  of  material  properties 
in  multimaterial  cells  is  given  in  Section  4.6.  The  following 
is  a summary. 

The  partial  volumes  of  the  materials  in  a multimaterial 
cell  are  determined  when  material  pressures  are  equilibrated. 

The  densities  of  the  materials  are  adjusted  in  an  iteration 

until  the  pressures  of  the  materials  equilibrate.  (See  Section  4.6.1) 

The  shear  yield  strength  of  a multimaterial  cell  is  a volume 
weighted  average  of  the  shear  yield  strengths  of  its  constituents. 

The  code  assumes,  however,  that  a multimaterial  cell  has  no 
tensile  strength,  i.e.,  negative  pressures  are  set  to  zero 
in  multimaterial  cells.  (See  Section  4.6.2) 

The  internal  energy  increment  of  a multimaterial  cell  is 
partitioned  by  partial  masses  in  SPHASE  and  HPHASE.  The 
internal  energy  is  updated  in  TPHASE  to  account  separately 
for  the  transport  of  internal  energy  and  the  thermalization 
of  kinetic  energy  of  each  material  in  the  cell.  (See  Section 
4.6.3) 

4.5.2  Free  Surface  Cells 

Any  cell  that  is  intersected  by  the  void  interface  is 
regarded  as  a free  surface  cell,  and  is  flagged  via  the  RHO 
array:  RHO  (NVOID,m)  = 1, where  NVOID  is  the  number  of  material 

* The  MFLAG  convention,  however,  imposes  a limit  of  99  material 
packages  in  a given  problem. 
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packages  plus  one  (NMAT+1),  and  m = MFLAG(k) -100.  A free 
surface  cell  can  contain  more  than  one  material  as  well  as 
the  void,  although  frequently  free  surface  cells  do  contain 
only  one  material. 

In  CDT,  if  the  equilibrated  pressure  of  a multimaterial, 
free  surface  cell  is  less  than  PMIN,  a user-controlled  in- 
put number  (See  Section  7.2.1),  the  densities  of  the  materials 
in  the  cell  are  not  modified  unless  the  sum  of  their  partial 
volumes  exceeds  the  volume  of  the  cell.  This  rule  prevents 
material  densities  in  a free  surface  cell  from  being  based 
on  the  mass  being  extended  over  the  entire  cell  volume.  How- 
ever, if  the  free  surface  cell  has  only  one  material,  the 
pressure  iteration  is  bypassed  and  the  density,  used  in  the 
equation  of  state,  RHOW,  is  the  total  cell  mass  divided  by  the 
total  cell  volume.  If  the  cell  pressure  based  on 
this  density  is  negative,  it  is  set  to  zero  and  the  density 
of  the  material,  RHO,  is  redefined  only  if  it  is  less  than 
the  assumed  density,  RHOW.  (See  also  Section  4.2.2.) 

In  SPHASE,  the  free  surface  cells  are  assumed  to  have 
no  strength.  Therefore  the  deviator  and  shear  stresses  of  a 
free  surface  cell  are  zero. 

In  DMCALC,  when  computing  the  mass  transported  from  a 
free  surface  cell  into  a non-free  surface  cell,  the  acceptor 
cell  density  instead  of  the  donor  cell  density  is  used,  since, 
in  general,  the  densities  of  the  non-free  surface  cells  are 
better  determined. 

4.5. ?>  Slipline  Cells 

Sliplines  in  HELP  are  not  restricted  to  grid  lines  or 
cell  diagonals.  They  do  coincide  with  material  interfaces, 
however.  Therefore  a cell  that  contains  a slipline  is  also 
an  interface  cell. 
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The  materials  in  a slipline  cell  are  accelerated  in 
inverse  proportion  to  their  densities  in  HPHASE.  (See 
Section  V.)  In  SPHASE  the  assumption  is  made  that  a slip- 
line cell  has  no  strengtli,  therefore  its  deviator  and  shear 
stresses  are  set  to  zero. 

The  velocity  components  of  the  "slave"  and  "master" 
materials  in  a slipline  cell  are  different.  The  cell 
velocity  components  represent  a mass  weighted  average  of 
the  material  velocity  components.  These  cell  velocities 
are  used  in  computing  the  work  terms  in  HPHASE  and  SPHASE 
and  in  moving  the  material  and  passive  tracers  in  MOVTCR. 

The  material  velocities  are  used  to  compute  the  mass  and 
momentum  transport  terms  in  TPHASE. 

4.6  DETERMINATION  OF  MATERIAL  PROPERTIES  IN  MULTIMATERIAL  CELLS 
Pressure  and  Density  Iteration  for  Multimaterial  Cells 

For  one-material  cells,  it  is  possible  to  compute  the 
pressure  directly  from  the  material  density  and  specific 
internal  energy.  However,  for  cells  containing  more  than  one  material, 
it  is  necessary  to  make  additional  assumptions  or  observations 
in  order  to  determine  the  cell  pressure.  In  the  present  method 
we  determine  the  cell  pressure  for  multimaterial  cells  by 
an  iteration  procedure.  Specifically,  the  densities  of  the 
various  materials  within  the  cell  are  varied,  subject  to  the 
condition  that  the  cell  be  exactly  filled  by  the  masses 
therein,  until  the  individual  pressures  converge  to  a common 
value,  which  is  taken  to  be  the  cell  pressure.  Specific  internal 
energies  are  held  constant  during  the  iteration.  This  process 
gives,  in  addition  to  the  cell  pressure,  the  densities  of  the 
individual  materials  for  subsequent  use  in  the  DMCALC  and 
TPHASE  transport  calculations. 

The  procedure  for  determining  the  pressure  and  the 
densities  in  a multimaterial  cell  is  actually  divided  into 
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two  parts:  a pre- iteration  calculation  to  insure  that  the 
cell  is  filled  exactly;  and  an  iteration  calculation  to  con- 
verge on  the  desired  P*s  and  p's. 


4. 6. 1.1  Pre- iteration  Calculation 

In  general,  the  masses  and  densities  which  exist  at  the 
beginning  of  the  calculation  do  not  exactly  fill  the  cell, 
since  cell  masses  are  changed  in  the  TPHASE  calculation  of 
the  previous  cycle.  The  pre- iteration  calculation  fills  the 
cell  exactly,  taking  into  account  the  (computed)  compressi- 
bilities of  the  various  materials  in  the  necessary  alteration 
of  densities.  Stepwise,  the  code 

1.  Computes  P^  from  the  old  value  of 
for  each  material  within  the  cell  (E^ 
is  assumed  constant  throughout) . 

2.  Varies  the  by  one  percent  to  compute 
new  P^. 

3.  Determines  the  constant-energy  compres- 

2 

sibilities,  called  for  each  material 

C?  = AP./Ap.  (4.10) 

where  the  AP^  and  Ap^  are  given  from 
the  results  of  (1)  and  (2). 

4.  Normalizes  (see  proof  below)  to  fill 
cell  exactly: 

AV^  = -AP/h?  (4.11) 


where 
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if  i 


i 
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t- 


It- 


= ( Pj  Ci  ) 


(4.12) 


and 


AP 


(VOL-VCELL) 


(4.13) 


VOL 


(4.14) 


with  "=  mass  of  constituent  i,  and 


VCELL  = volume  of  the  cell. 


With  the  above  increments  in  specific 
volume,  the  new  specific  volumes  are 
given  by 


" V.  , + AV. 
1 i,old  1 


(4.15) 


Then, 


Pi  . l/Vi 


(4.16) 


gives  the  desired  nev;  densities  which 
cause  the  cell  to  be  filled  exactly. 


It  is  easy  to  see  that  the  calculation  defined  above 
causes  the  cell  to  be  filled  exactly.  To  show  this. 
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Summing , 


’’i/Pi.old  - "'i 

m. (VOL-VCELL) 

m . / D • n J “ — 
i'  ,old 


= VOL 


■E 


m , 


>'1 


- (VOL-VCELL) 


(4.17) 


= VCELL 


(4.18) 


showing  tliat  the  new  volum^js  m^V^  of  the  constituents  add  up 
to  the  total  cell  volume,  as  desired. 


4 . 6 . 1 . 2 Iteration  Calculation 

Given  P.  = £■ (v. , E.)  for  materials  i=  1,...,M,  the  cell 
1 1 \ 1 1/ 

pressure  P = P-  = ?2  = ...  = Pj^  is  found  by  varying  the  V^ 
until  the  P^  are  equal  within  some  specified  accuracy, 
subject  to  the  conditions  that  the  cell  remain  exactly 
filled  and  that  the  are  constants. 

The  equation  for  volume  conservation  is 


J^m^  AV^  = 0 (4.19) 

and  the  equation  that  causes  material  i to  undergo  a change 
in  specific  volume  AV^ , such  that  its  pressure  changes  from 
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its  current  value  = ^i(^i’^i)  ^ value  common 

to  all  materials  in  the  cell  at  the  end  of  the  iteration 
cycle,  is 


(4.20) 


These  equations,  for  i = 1,...,M  can  be  regarded  as  M+1 
equations  for  the  M+1  quantities  AV^  and  F"  . The  other 
quantities  in  the  equations  are  either  known  constants  (m^) 
or  are  updated  each  cycle  of  the  iteration  [the  and  the 
|3V^/9P^j  ] and  are  taken  to  be  constants  while  solving 

^i 

— ^ . *1 

these  equations  for  AV^  and  P . The  solutions  are 


where 


rn+1 


AV. 

1 


p.  w” 

LmJ  1 1 

(p?  - 1 


1 

FT 

1 


rn+1 
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(4.22) 


(4.23) 


w.  = m./h. 
1 1 1 


(4.24) 


That  these  equations  are  solutions  to  the  given  equations  can 
be  verified  by  direct  substitution. 

The  iteration  would  be  exact  (single  cycle  convergence) 
if  the  input  coefficients  (c'P./9V.)  were  constants,  since 


the  solution  does  not  involve  additonal  approximations.  In 
the  iteration,  these  coefficients  are  refined  each  step  of 
the  iteration  by  making  use  of  the  P^,  points  which  were 
generated  during  the  previous  step.  Further,  some  pains  are 
taken  to^  determine  realistic  starting  values  of  these 
coefficients,  in  the  pre- iteration  calculation,  in  the  belief 
that  this  procedure  will  hasten  convergence  sufficiently  to 
reduce  total  computing  time. 

Stepwise,  the  iteration  calculation  is 

as  follows: 

1.  Compute  new  and  AV^  by  the  above 

formulas,  using  known  values  of  and 

(av,/aP,)^^. 

n 4»  1 

2.  Compute  new  v”  by  adding  the  above  AV^ 
to  the  old  specific  volumes. 

3.  Compute  new  P^  , using  the  updated  V^. 

4.  Compute  refined  values  of  the  coefficients 
|3V^/3P^j  by  using  this  last  point  [from 

(2)  and  (3)  above]  and  that  previous 
point  (two  such  points  are  saved  for  each 
material  during  the  iteration)  which  is 


The  iteration  is  completed  when  the  all  equal  P to  some 
preset  accuracy,  PRCNT  (See  Section  7.2.1).  In  most  exper- 
ience to  date,  convergence  has  been  obtained  in  two  or 
three  steps,  where  a convergence  criterion  |P^  - F|  < 10  P 
was  used. 

4. 6. 1.3  Equation  o£  State  Modifications  for  the  Iteration 


roceoure 

The  Tillotson  equation  of  state  (described  in  Section 
2. 3. 1.2)  has  been  used  in  most  work  to  date  with  HELP  and  is 
sufficiently  general  to  provide  a satisfactory  description  of 
most  media  which  is  initially  condensed  and  which  may  vaporize 
as  a result  of  heating.  However,  for  some  materials  which  are 
sufficiently  expanded  and  cold  (E  < E^),  the  slope  of  a con- 
stant energy  curve  (in  the  P-V  plane)  is  positive.  For  pure 

cells,  these  meaningless  states  are  avoided  by  replacing  the 

2 

quantity  (Ap  + Bp  ) with  its  minimum  value  (a  negative  pressure) 
if  the  cell  density  is  less  than  the  density  corresponding  to 
this  minimum  value. 

To  assure  convergence  of  the  pressure  iteration  in 
multimaterial  cells,  it  is  necessary  to  have  a continuous 
equation  of  state  and  for  8P/8V  to  be  everywhere  negative. 

To  this  end  the  term 


Ap  + Bp‘ 


(4.26) 


is  replaced  by 


Ap  + h(v)Bp‘ 


where 


1.  h(v)  ■ +1  for  V £ Vjj  (compressed 
states) 

2.  h(v)  decreases  linearly  to  -1  in 
the  interval 


^0  < V < V^  ^ 


(4.28) 


3.  h(v)  - -1  for  V > V 


- ’'o  ina 


The  effect  of  these  modifications  is  to  change  B to  -B  in  a 
continuous  manner  and  hence  insure  that  8P/3V  is  always  nega- 
tive. Negative  pressures  that  may  result  after  the  iteration 
is  completed  are  set  to  zero,  since  tensile  states  in  inter- 
face cells  are  assumed  not  to  exist.  In  view  of  this  latter 
assumption,  the  indicated  equation  of  state  changes  do  not 
normally  affect  the  physical  result,  since  they  alter  only 
states  that  would  ordinarily  be  tensile.  The  purpose  of  the 
changes  is  to  make  the  equation  of  state  well-behaved,  so 
that  the  iteration  converges  properly. 


4.6.2  Shear  Yield  Strength  of  Multimaterial  Cells 


The  shear  yield  strength  of  a multimaterial  cell  is  a 
volume  weighted  average,  7,  of  the  shear  yield  strengths- of 


its  constituents,  Y„: 

n 


(4.29) 


If,  however,  one  of  the  constituents  is  an  ideal  gas  or  a 
high  explosive,  (MAT(n)  ^ 20),  the  shear  yield  strength  of 


2.. 
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the  cell  is  set  to  zero.  Likewise,  if  the  cell  contains  the 
void  interface  (RHO(NVOID,m)  ■ 1)  or  a slipline  interface 
(THETA (m)  >0),  with  m ■ MFLAG(k)  - 100,  the  shear  yield 
strength  of  the  cell  is  set  to  zero.  And  finally,  the  shear 
yield  strength  of  the  cell  is  zero  if  the  density  of  any 
one  of  its  constituents  is  below  the  density  corresponding 
to  the  maximum  tensile  strength  of  the  constituent 


(RHO(n,m)  <p  *AMDM(n)).  This  presumes  that  a material  that 


n 


has  failed  in  tension  has  no  shear  yield  strength. 


4.6.3  Partitioning  of  Internal  Energy  Increments  in 
Multimaterial  Cells 


The  internal  energy  increments,  AEj,  computed  by  HPHASF 
and  by  SPHASE  are  partitioned  amongst  the  materials  in  a 
multimaterial  cell  in  proportion  to  their  masses.  For  ex- 
\ ample,  the  internal  energy  increment  for  material  n is 


AEr  »AE 


'I 


n 


I 


■vS 


rtii 


(4.30) 


Iii^  other  words  the  HPHASE  arid  SPHASE  gain  in  specific  internal 
energy  is  the  same  for  all  the  materials  in  the  cell. 

The  situation  is  somewhat  more  complex  in  TPHASE,  since 
the  iftternal  energy  increment  is  due  not  only  to  transport 
but  also  to  the  thermal ization  of  kinetic  energy.  It  is 
natural  to  imagine  that  the  interacting  materials  (all  of  the 
masses  of  material  n within  the  cell  at  the  end  of  the  time 
step)  have  collided  in  a coordinate  system  in  which  their 
center  of  mass  has  zero  velocity. 

First,  we  recall  that  the  new  velocity  components  for  , 
material  n are  chbsen  to  conserve  momentum 
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and  that  the  new  mass  of  material  n is 


■ bI|P  + iC  ^™4n 

n n j"  in 


(4.33) 


Here  the  subscript  i denotes  the  four  sides  of  the  cell;  the 
subscript  n denotes  the  material  package,  and  the  superscript 
(n)  denotes  pre-TPHASE  quantities.  The  velocity  associated 
with  each  of  the  transport  masses  is  the  donor  cell  material 
velocity;  e.g.,  if  a mass  Am^^  entered  the  cell  under  consider- 
ation through  side  1,  the  associated  velocity  u^^^  is  the 
velocity  of  material  n in  the  cell  from  whence  the  mass  came. 

We  note  first  that  as  determined  to  con- 

n ' n 

serve  momentum  by  Eqs.  (4.31)  and  (4.32)  above,  are  the 
velocities  of  the  center  of  mass  of  the  material  n in  the  cell 
at  the  end  of  the  time  step.  This  means  that  a component  mass 
with  pre-TPHASE  velocities,  following 

amount  of  kinetic  energy  thermalized 


,(n+l)  . ^(n) 


in  the  center  of  mass  coordinate  system.  Thus  the  total  ther- 
malized KE  for  material  n is 
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,(n+l) 


(n+1) 


Some  interpretation  of  this  equation  clarifies  its  meaning: 
the  first  term,  plus  those  terms  from  the  summation  for  which 
6min  < 0 (material  leaving  the  cell  with  donor 

cell  velocities  represents  the  ther- 

in  n ’ in  n 

malized  KE  for  that  fraction  of  the  material  which  remains  in 

the  cell.  The  additional  terms  (6m.  > 0 in  the  summation) 

in 

represent  the  thermalized  KE  associated  with  that  material 

which  came  into  the  cell  during  the  time  step.  For  these 

terms  the  uf^^  v|”^  refer  to  the  velocities  of  material  in 
in  in 

the  neighbor  cells  which  acted  as  donors. 

The  change  in  internal  energies  is  again  given  by 


EaEii„  . TKE„ 


(4.35) 


In  other  words,  material  n gains  an  amount  of  total  internal 
energy  which  is  that  due  to  transport,  plus  the  amount  of 
thermalized  kinetic  energy  for  material  n. 

The  above  discussion  assumes  that  the  material  velocities 
are  not  equal  to  the  cell  velocities,  i . e. , that  the  cell  con- 
tains the  slipline.  If  the  cell  does  not  contain  a slipline, 
equations  (4.33)  through  (4.35)  still  hold,  only  the  material 
velocities  are  defined  using  cell  quantities  as  follows: 
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CHAPTER  V 


SLIPLINES 


S . 1 GENERAL  COMMENTS 


Sliplines  have  existed  in  Lagrangian  codes  for  several 
years  ani  some  of  our  thinking  and  terminology  comes  from 
these  earlier  accomplishments.  For  example,  in  HELP  one 
package,  is  considered  the  "master"  along  which  another  "slave" 
package  slips.  However,  an  important  difference  between  this 
current  effort  using  an  liulerian  code  and  the  earlier  Lagrangian 
models  is  that  the  sluipe  of  the  slip  surface  is  determined  by 
th^  stress  field  and  need  not  coincide  with  grid  lines  or  node 
points. 

The  slipline  in  HELP  is  also  a material  package  inter- 
face; therefore,  the  slipline  cells  arc  also  interface  cells. 

The  angle  of  the  slipline  in  a given  slipline  cell  is  measured 
from  the  x-axis  of  the  grid,  using  *the  intercepts  of  the  slip- 
line interface  with  the  cell  boundaries.  Tiiese  intercepts 
are  computed  while  processing  the  interfaces  in  TRACS  and  are 
stored  by  PTSAV.  The  stored  intercepts  are  then  used  by 
THETAS  to  compute  the  angle  of  the  slipline  in  each  slipline 
cell. 


Since  the  slave  and  master  materials  in  a slip  cell 
develop  different  velocity  components  tangent  to  the  slipline, 
it  is  necessary  to  store  distinct  velocity  components  for 
each  material.  (In  an  earlier  version  of  HELP,  it  was  assumed 
that  all  materials  in  an  interface  celJ  would  have  the  same 
velocity  components.)  Two  material  arrays,  US,  \^S,  have  been 
added  to  store  the  material  velocities  in  interface  cells. 

To  describe  the  angle  of  the  slipline  across  a given  slip 
cell,  the  THETA  array  was  added.  When  THETA (m)=-l  (m  = MFLAG(k) 
100),  the  interface  cell  does  not  contain  the  slipline;  when 
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I T[!nTA(m)  > 0 the  slipline  docs  cross  thr-  interface  cell  at 

i an  anslc  (in  radians)  equal  to  THETA (m). 

I It  should  be  noted  that  more  than  one  material  package 

I • can  be  designated  as  a slave  or  master.  Therefore  it  is 

I possible  for  a sliplinc  cell  to  contain  more  than  one  slave 

i ~ or  master  package.  In  that  case,  the  velocities  of  all  the 

• slave  materials  in  the  cell  arc  equilibrated  as  are  the 

I 

velocities  of  all  the  master  materials.  The  slave  density, 

, Pg,  used  in  the  equations  that  follow  is  an  average  of  the 

densities  of  all  the  slave  materials  in  a given  slip  cell, 
likewise  for  the  master  density, 


5.2  RELATIVE  ACCELERATION  OF  MASTER  AND  SLAVE 

MATERIALS  rUVCALO 

A cell  that  contains  the  sliplinc  has,  in  a sense,  failed. 
Therefore  the  code  assumes  the  material  in  a slip  cell  has 
no  shear  yield  strength  and  the  deviator  stresses  of  a slip 
cell  are  set  to  zero.  The  only  stresses  acting  on  a slip 
cell  are  simply  the  hydrostatic  pressures.  The  acceleration 
of  slip  cells  due  to  stresses  is  computed  entirely  in  HPHASE. 
After  updating  the  momenta  of  a slipline  cell,  IIPHASE  calls 
UVCALC  to  establish  the  relative  acceleration  of  the  master 
and  slave  materials.  Given  the  post-MPMASE  cell  centered 
momenta  (m(k)u(k),  m(k)v(k))  and  the  angle  9 of  the  slipline 
through  cell  k,  the  master  and  slave  material  velocities,  uSj„» 
vSjj^,  and  us^,  vs^.,  are  updated  by  solving  the  following  four 
simultaneous  linear  equations: 

vs^  + m„vs  = (m^  + m„)  v(k)  (5.1) 

ssmm  sm"^^' 

ni..  + m„us„  = (m^  + m„)  u(k)  (5.2) 

ssmm^sm'^^'^  ^ ' 
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s in  0 
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uSg  sin  0 - vSg  cos  0 = u 


% 


vs  cos  0 
m 


Ps  [(uSy-us^)  cos  0 + (vs^  - vs^)  sin  0] 


(5.4) 


where  usl , vs'.,  us',  vs'  arc  pre-llPHASl:  slave  and  master 
material  velocities.  The  updated  material  velocities  must 
conserve  the  cell's  axial  and  radial  momenta  (nqs.  (5.1)  and 
(5.2)).  The  m.aster  and  slave  velocity  components  normal  to 
the  slipline  must  be  eciiial  fP.q.  (5.3)),  and  the  HI’HASE  change 
to  the  slave  and  master  velocity  components  tangent  to  the 
slipline  must  he  inversely  proportional  to  the  slave  and 
master  material  densities,  p^.,  (Eq.  (5.4)). 

Solving  these  four  equations  results  in  the  thermalization 
of  some  of  the  slip  cell's  kinetic  energy.  This  thermalized 
kinetic  energy  is  converted  to  internal  energy  and  is  partitioned 
amongst  the  materials  in  the  cell  in  proportion  to  their  masses. 


5.3  TRANSPORT  OF  MASTER  AND  SLAVE  NIATERIALS 

The  mass,  momentum  and^  energy  transport  terms  are  com- 
puted in  DMCALC  and  TPHASE.  Since  some  interface  cells, 
namely  those  along  the  slipline,  have  different  velocity 
components  for  the  master  and  slave  materials,  the  mass  and 
momentum  transport  equations  for  interface  cells  must  employ 
material  velocities  rather  than  cell  velocities. 

For  example,  the  mass  of  material  n to  be  transported 
across  the  i cell  face  is 


6mn 


At 


(5.5) 
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where  p^  is  the  density  of  material  n in  the  donor  cell,  u5 
is  the  mass  transport  velocity  based  on  the  material  n veloc- 
ities of  the  donor  and  acceptor  cells,  and  is  the  i cell 

^i 

face  area  associated  with  material  n.  Likewise,  the  momentum 
transport  terms  are  based  on  material  velocities  rather  than 
cell-centered  velocities: 

«mun  “ % * 

^mVn  - 6m^.  vs^  . C5.7) 

Since  the  material  velocities  equal  the  cell  velocities 
in  cells  that  do  not  contain  a slipline,  the  equations  a re 
result  in  the  same  rnasc  and  momentum  transport  for  non- slip 
interface  cells  as  do  equations  that  use  cell-centered  quan- 
tities . 


5.4  MODIFICATION  OF  POST-TPHASE  M^VSTER  AND 

SLAVE  VELOCITIES  rUVMODl 

After  the  material  velocities  of  a slip  cell  have  been 
updated  in  TPHASE,  they  must  be  altered  to  insure  that  the 
material  velocity  components  normal  to  the  slipline  are  equal. 

Subroutine  UVMOD  is  called  from  TPHASE  after  all  cells 
have  been  updated  for  the  effects  of  transport.  UVMOD  searches 
for  interface  cells  that  contain  a slipline  (THETA(m)  > 0)  and 
modifies  the  post-TPHASE  master  and  slave  material  velocities 
of  those  cells  so  that  their  components  normal  to  the  slipline 
are  equal,  i.e.,  so  that  Eq.  (5.3)  in  the  discussion  above  is 
satisfied.  The  material  velocities  of  non-slip  interface  cells 
are  not  modified  by  UVMOD. 
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5.5  USE  OF  SLIPLINES 


Sliplines  in  HELP  have  been  employed  in  a variety  of 
situations;  for  example,  at  the  high  explosive-metal  inter- 
faces of  shaped  charge  and  fragmenting  munitions  calculations, 
and  at  the  interface  between  the  plug  and  target  when  modeling 
thin  plate  perforation. 

When  there  is  a large  density  difference  between  the 
master  and  slave  materials  (as  between  an  expanded  HE  and  a 
metal  casing) , the  relative  motion  of  the  two  packages  is 
generated  primarily  by  the  HPHASE  accelorations> which  are 
density  weighted  and  which  satisfy  Eqs.  (5.1)  through  (5.4) 
in  the  discussion  above.  However,  when  the  densities  of  the 
master  and  slave  materials  are  close  (as  between  a plug  and 
a target),  the  relative  motion  of  the  packages  results  primarily 
from  the  different  momentum  transported  to  the  master  and 
slave  materials  in  TPHASE.  For  example,  given  that  the  pro- 
jectile and  plug  are  both  slave  packages,  and  the  target  is 
a master  package,  the  material  velocities  of  the  plug  and 
projectile  in  a slip  cell  are  equilibrated,  thereby  trans- 
ferring the  projectile  momentum  to  the  plug  but  not  to  the 
target. 
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CflAPTER  VI 
PLUGGING  FAILURE 


A plugging  model  has  been  incorporated  into  the  HELP 
code  which  the  user  activates  if  he  or  she  anticipates  plug- 
ging failure.  As  the  following  discussion  will  clarify,  the 
decision  to  use  the  plugging  failure  model  must  be  made  when 
the  calculation  is  generated.  As  the  model  now  exists,  its 
application  is  restricted  to  the  impact  of  right  circular 
cylinders  or  truncated  cones  into  homogeneous  targets  of 
finite  thickness.  Furthermore,  the  user  must  supply  an 
empirically  derived  material  parameter  (PLWMIN)  which  deter- 
mines the  shape  and  mass  of  the  plug  that  is  formed 

In  the  sections  that  follow,  the  procedures  for  generating 
a plugging  calculation  as  well  as  the  code’s  method  for  evolving 
the  plug  will  be  discussed. 

6.1  GENERATION  OF  A PLUGGING  CALCULATION 

6.1.1  Designation  of  Material  Packages 

The  plugging  model  assumes  the  projectile  is  package  1, 
the  plug  is  package  2,  and  the  target  is  package  3.  The  pro- 
jectile can  be  made  up  of  several  m.itcTial  packages  as  long 
as  the  part  of  the  projectile  which  will  interact  with  the 
target  is  homogeneous  and  is  designated  as  package  1.  Further- 
more, the  projectile  and  plug  packages  must  be  "slave"  packages 
and  the  target  a "master"  package. 


6.1.2  Placement  of  Plug  Tracer  Particles 


Unlike  the  other  material  packages  in  a HELP  calculation, 
the  dimensions  of  the  plug  package  are  unknown  at  the  beginning 


of  the  calculation.  However,  the  user  does  specify  in  the 
input  deck  a set  of  plug  tracer  particles  which  later  in  the 
calculation  are  automatically  redefined  so  as  to  circumscribe 
the  plug  package  as  it  evolves.  Initially  these  plug  tracers 
coincide  with  the  projectile  and  target  tracers  along  the 
projectile-target  interface.  These  initial  tracer  positions 
are  based  on  the  assumption  that  the  vertical  edge  of  the 
plug  package  begins  at  the  top  right  edge  of  the  projectile. 

In  the  example  below,  the  relationship  between  the  pro- 
jectile, target  and  plug  tracer  particles  at  time  T = 0 is 
illustrated : 


Projectile  - The  last  six  projectile  tracer 
particles  (excluding  the  dummy  end  point*) 
are  at  positions  F through  A. 


* The  last  tracer  of  a package  is  a dummy  point.  It  is  used 
only  to  indicate  the  package  has  been  completely  circum- 
scribed. 
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Target  - The  first  six  target  tracer  particles 
are  at  positions  A through  F. 

Plug  - The  first  six  plug  tracer  particles  are 
at  positions  A through  F.  The  final  six  plug 
tracer  particles  (excluding  the  dummy  point) 
are  at  positions  F through  A. 


w. 


There  is  a set  of  n duplicate  plug  and  target  tracer  particles 
placed  at  position  F in  order  to  define  the  vertical  edge  of  the 
plug  package  as  it  is  extended.  (See  Section  6.2. 2.1.)  These 
n duplicate  particles  are  generated  automatically  by  PLGADD, 
given  the  user- supplied  definition  of  the  plug  slipline  end- 
points (nbGSI)(2),  NENDSD(2)j.  The  number  of  duplicate 
particles  generated  is  (NENDSD(2)  - NBGSD(2)  + 1)  and  is 
therefore  controlled  by  the  user. 

6.1.3  Placement  of  Passive  Tracer  Particles 

The  user  must  specify  a region  of  the  target  in  which  the 
plug  is  expected  to  form.  In  each  cell  of  that  region,  four 
passive  tracer  particles  are  automatically  placed  by  SETUPA. 

These  passive  particles  are  used  to  track  material  within  the 
target  , to  accumulate  its  plastic  work,  and  to  compute  its 
angle  of  maximum  shearing  stress.  This  information  is  needed 
for  the  plugging  failure  criterion  discussed  in  Section  6.2.1. 

6.1.4  Slipline  Specification 

To  simulate  the  action  of  a plug  being  pushed  out  of  a 
target  plate,  it  is  necessary  to  generate  slip  surfaces  be- 
tween the  plug  and  the  target  and  between  the  deforming  projectile 


head  and  target  lip.  These  slip  surfaces  do  not  exist,  of 
course,  until  the  plug  evolves.  In  a plugging  calculation 
the  tracer  particles  that  define  the  endpoints  of  the  slip- 
line for  each  package  are  all  initially  at  the  top  right 
corner  of  the  projectile,  where  the  plug  begins  to  form. 

6.2  METHOD  OF  EVOLVING  THE  PLUG 

To  simulate  plugging  failure,  the  HELP  code  essentially 
evolves  a plug  package  the  shape  and  mass  of  which  are  deter- 
mined by  a plastic  work  failure  criterion  and  by  the  assoc- 
iated angle  of  maximum  shearing  stress. 

6.2.1  Plugging  Failure  Criterion 

The  failure  criterion  used  to  propagate  the  vertical 
edge  of  the  plug  is  based  on  the  specific  plastic  work  of  the 
target  material  in  the  region  near  the  endpoint  of  the  plug 
edge.  The  plastic  work  of  this  material  is  approximated  by 
using  passive  tracer  points  that  move  with  the  target  material 
and  by  incrementing,  by  an  area  weighted  average,  the  plastic 
work  of  the  material  associated  with  each  point.  Also  assoc- 
iated with  each  point  on  a given  cycle  is  an  area  weighted 
average  angle  of  maximum  shearing  stress.  Once  the  plastic 
work  associated  with  one  or  more  points  in  the  region  of  the 
endpoint  exceeds  a predetermined  cutoff  value  (PLWMIN) , the 
average,  current  angle  of  maximum  shearing  stress  associated 
with  those  points  is  used  to  determine  the  direction  in  which 
the  vertical  edge  of  the  plug  is  extended. 

Each  cycle  in  SPHASE,  a plastic  work  increment  and  an 
angle  of  maximum  shearing  stress  are  associated  with  the 
center  of  cells  in  the  target.  An  area  weighted  average  of 
these  cell-centered  quantities  is  used  to  define  an  angle  of 
maximum  shearing  stress  and  to  update  the  total  plastic  work 
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for  the  material  associated  with  each  moving  point.  For 
example,  consider  point  j in  the  figure  below. 
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The  plastic  work,  W,  associated  with  point  j at  time  n is 
defined  as 


w-  - 


(6.1) 


The  average  angle  of  maximum  shearing  stress,  a,  for  the 
point  is  defined  as: 


“j  ’ • 


(6.2) 


1 


In  both  cases  A.  denotes  the  overlap  area  of  the  pseudo-cell, 
1 11  n • 

centered  at  point  j,  with  cell  i.  AW^  and  are  the  time  n 

plastic  work  increment  and  angle  of  maximum  shearing  stress, 

respectively,  associated  with  the  center  of  cell  i. 

Each  cycle,  after  the  plastic  work  associated  with  the 
moving  points  has  been  updated,  the  code  considers  a region 


which,  as  illustrated  below,  is  one  cell  in  area  and  lies  above 
the  intercept  between  the  topmost  horizontal  grid  line  and  the 
plug's  vertical  surface.  The  radial  position  of  the  region  is 
centered  about  this  intercept. 


If  one  or  more  of  the  moving  points  which  lie  within  that 
region  has  an  accumulated  plastic  work  total  greatei  than  the 
specified  empirical  value,  PLWMIN,  then  the  vertical  edge  of 
the  plug  is  extended  to  the  next  horizontal  grid  line.  The 
angle  at  which  the  plug  edge  is  extended  is  an  average  of  the 
current  angles  of  maximum  shearing  stress  associated  with 
those  points  within  the  specified  region  that  have  exceeded  the 
plastic  work  cutoff. 

For  example,  in  the  case  illustrated  above,  the  code  would 
consider  the  plastic  work  associated  with  passive  tracer  particles 
c,  d,  and  f,  only,  since  passive  tracers  a,  b and  e are  outside 
the  specified  region.  If,  on  a given  cycle,  the  plastic  work 
associated  with  all  three  points  c,  d and  f exceeds  the  specified 
value,  PLWMIN,  the  angle, a,  at  which  the  plug  will  be  extended 
will  be  a = ^cl  ^ “f^  ' point  c satisfied  the 

failure  criterion,  then  u = a^. 
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When  the  criterion  for  extending  the  plug  surface  is 
met,  the  plug  package  itself  must  be  expanded.  The  mechanisms 
for  this  expansion  are  described  in  the  following  sections. 

6.2.2  Conventions  Governing  Plug  Tracer  Particles 

There  are  three  stages  of  the  plug  formation.  The  con- 
ventions governing  the  initial  placement  of  the  plug  tracer 
particles  have  been  given  in  Section  6.2.1.  The  conventions 
for  the  intermediate  and  final  stages  are  given  below. 


6. 2. 2.1  Partially  Formed  Plug 

Once  the  plug  edge  is  extended,  the  projectile,  target, 
and  plug  tracer  particles  follow  the  conventions  illustrated 
by  the  example  below,  in  which  5 of  the  set  of  n duplicate  plug 
package  tracer  particles  generated  initially  by  PLGADD  (see 
Section  6.1.2)  now  define  the  vertical  edge  of  the  plug. 


I 


Proj ectile  - The  last  six  projectile  tracers 
(excluding  the  dummy  tracer)  are  at  positions 
F through  A. 

Target  - The  first  six  target  tracers  are  at  P 
through  K.  The  next  (n-5)  are  also  at  K,  where  n 
is  the  number  of  duplicate  points  initially  speci- 
fied by  NBGSD(2)  and  NENDSD(2).  The  next  five 
are  at  J through  I'. 

Plug  - The  first  six  plug  tracers  are  at  A through 
F.  The  next  four  are  at  G through  J.  The  next 
(n-5)  are  at  K,  and  the  last  five  (excluding  the 
dummy  tracer)  are  at  L through  P. 

Each  time  the  plug  is  extended  all  but  one  of  the  dupli- 
cate particles  is  moved  to  the  new  endpoint  of  the  vertical 
edge  of  the  plug.  The  one  particle  left  behind  helps  to  de- 
fine the  vertical  plug  interface. 

It  should  be  noted  here  that  the  interface  in  the 
example  above  does  not  represent  a discontinuity  in  the 
material  as  does  the  slipline  interface  KF.  Since  FF  is  a 
psuedo- interface,  the  material  flow  field  is  unaffected  by 
its  presence. 


6. 2. 2. 2 Completely  Formed  Plug 

Once  the  plug  edge  is  extended  to  the  free  surface  back 
of  the  target  and  the  projectile  and  target  lips  coincide,  the 
projectile,  target  and  plug  tracer  particles  follow  the  con- 
ventions outlined  in  the  example  below. 
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Projectile  - The  last  six  projectile  tracers 
(excluding  the  dummy  tracer)  are  at  positions 
F through  A.  The  target  and  projectile  have 
joined  between  F and  F'  and  the  slipline  extends 
from  P to  F'. 

Target  - The  target  package  no  longer  connects 
with  the  axis.  The  first  eleven  tracers  are 
at  P through  F.  The  last  target  tracer  (excluding 
the  dummy  tracer)  is  also  at  P in  order  to  completely 
circumscribe  the  target  package. 

Plug  - The  first  six  plug  tracers  are  at  A 
through  F.  The  next  ten  are  at  G through 
P,  and  the  last  four  (excluding  the  dummy 
tracer)  are  at  P through  S.  Note  that  there 
are  two  duplicate  plug  points  at  P;  one  is 
moved  with  plug  velocities  and  defines  the 


top  right  corner  of  the  plug  as  it  is  pushed  out; 
the  other  is  moved  with  target  velocities  and  defines 
the  top  left  corner  of  the  target.  The  last  three 
tracers  of  the  plug  are  coincident  with  what  were 
the  last  three  tracers  of  the  target  just  before  the 
plug  was  completely  formed. 

6.2.3  Movement  of  Plug  Tracer  Particles 

As  described  in  Section  4.1.2,  tracer  particles  are  moved 
with  a density  and  area-weighted  average  velocity.  However, 
to  better  simulate  plugging  and  penetration  mechanisms  this 
method  has  been  specialized  for  those  plug,  projectile  and 
target  tracer  particles  that  define  the  plug  package  once  it 
has  begun  to  evolve. 

1.  The  tracer  particles  along  the  projectile- 
plug  interface  are  not  moved  with  target 
material  velocities,  only  plug  and  pro- 
jectile velocities. 

2.  The  tracers  along  the  slipline  edge  of  the 
plug,  except  for  the  top  endpoint,  are  also 
moved  with  only  plug  and  projectile  veloci- 
ties. 

3.  The  tracers  along  the  top  surface  of  the  plug, 
except  the  top  endpoint,  are  moved  normally 
using  all  velocities  in  the  overlap  region. 

4.  The  top  endpoint  of  the  plug  surface  is  moved 
with  only  target  velocities. 

The  following  example  illustrates  the  way  in  which  a 
tracer  in  category  1 above  is  moved.  If  the  area  associated 


(2)  \lARdliT  1.1) 


1^1 

1 


kv\  kd 


ppolncTili;  (1 


with  one  of  these  tracers  overlaps  a cell  which  contains  only 
target  material,  e.g.,  cell  kd  in  the  figure,  the  overlapped  area 
and  its  associated  velocity  component  are  ignored  in  the  de- 
finition of  the  tracer's  velocity  components.  The  axial 
velocity  of  the  circled  point,  therefore,  is 


% ' (^k'^VSi  ♦ vS2)/2  * 


. A^|,-(vsj  . VSj)/2)/(A^  * A|^„  * A^^^) 


(6.3) 
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where  is  the  overlap  area  between  cell  i and  the  pseudo- 
cell, and  vs^  is  the  axial  velocity  of  material  j in  cell  i. 

Notice  that  the  target  velocity,  vs^,  in  the  cells  was  not  in- 
cluded in  the  average.  The  cell  centered  velocities  are 
used  in  the  averages  only  when  the  overlap  cell  is  pure.  The 
reasoning  behind  these  special  procedures  is  that  the  slipline  rep- 
resents a discontinuity  in  the  velocity  field  and  a tracer  on 
one  side  of  this  line  should  not  be  affected  by  the  motion  of 
the  material  on  the  other  side. 

The  tracers  along  the  free  surface  of  the  target  are 
moved  with  target  velocities  only.  Tracers  along  the  free 
surface  of  the  plug  are  moved  with  plug  velocities  only.  The 
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tracers  along  the  free  surface  of  the  projectile  are  moved 
with  both  plug  and  projectile  velocities  in  the  case  where 
the  pseudo  cell  overlaps  an  interface  cell  which  contains 
plug  material.  Otherwise  they  are  moved  only  with  projectile 
material  velocities. 

Once  the  plug  surface  has  extended  to  the  rear  free 
surface  of  the  target,  the  top  endpoint  of  the  slipline  is 
considered  a part  of  the  target  free  surface;  therefore,  it 
is  moved  with  target  material  velocities  only.  As  the  plug 
is  pushed  out,  the  slipline  tracers  below  the  top  endpoint 
will,  in  general,  have  a larger  axial  component  of  velocity 
than  the  endpoint  since  they  are  moved  with  plug  and  pro- 
jectile material  velocities.  The  plug  slipline  tracers  are 
accumulated  at  the  endpoint.  The  corresponding  target  slip- 
line tracers  are  removed  and  the  number  of  target  tracers  is 
reduced  accordingly. 

6.2.4  Conversion  of  Target  Material  into  Plug  Material 
6. 2.4.1  Redefinition  of  Tracer  Particles 

When  the  plug  surface  is  extended,  the  top  edge  of  the 
plug  is  redefined,  and  the  plug  package  boundary  encompasses 
a larger  volume  as  illustrated  below. 
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The  tracers  at  the  top  of  the  plug  are  moved  a fraction 
above  the  next  horizontal  grid  line,  which  is  also  the  new 
y-cnordinate  of  the  slipline  endpoint.  However,  as  illus- 
trated below,  those  tracers  at  the  top  edge  of  the  plug 
that  arc  already  above  the  new  slipline  endpoint  are  not 
moved. 
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Until  the  plug  is  completely  formed,  the  upper 
endpoint  of  the  slipline  is  on  the  top  horizontal  grid  line 
of  the  last  failed  row.  However,  the  top  and  vertical  edges  of  the 
plug  may  extend  above  this  horizontal  line,  as  shown  below, 
where  • denotes  the  slipline  endpoints. 
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However,  when  the  plug  edge  is  extended  into  a cell  that  con- 
tains the  free  surface,  the  rear  of  the  target  has  been 
reached,  and  the  plug  edge  is  extended  beyond  the  grid  line 
until  it  intersects  the  free  surface,  as  illustrated  below. 


The  tracers  at  the  top  of  the  plug  are  not  moved  in  this 
case.  Instead,  the  target  tracers  to  the  left  of  the  new 
plug  surface  endpoint  replace  the  tracers  at  the  top  edge 
of  the  plug.  At  this  time  in  the  calculation  the  number  of 
plug  and  target  tracers  is  reduced.  All  but  two  of  the  re- 
peated plug  tracers  at  the  plug  surface  endpoint  are  removed. 
One  of  these  will  define  the  top  right  corner  of  the  plug  as 
it  is  pushed  out;  the  other  will  define  the  top  left  corner 
of  the  target  and  will  also  be  the  endpoint  of  the  slipline. 

The  target  package  is  no  longer  attached  to  the  axis,  and  the 
first  and  last  target  tracers  define  the  new  slipline  endpoint. 

Accordingly,  two  free  surface  tracers  are  added  at  the 
slipline  endpoint;  one  will  move  with  the  plug,  the  other 
will  move  with  the  target  and  remain  at  the  slipline  endpoint. 
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6 . 2 . 4 . 2 Redefinition  of  Target  and  Plug  Masses 

When  the  plug  edge  is  extended  and  the  plug  package 
boundary  is  enlarged,  the  code  converts  some  of  the  target 
material  into  plug  material.  To  do  this,  the  code  computes 
the  volume  to  be  associated  with  plug  material  in  cells 
containing  the  plug  edge  extension.  Using  the  intercepts  of 
the  plug  edge  extension  with  the  cell  boundaries,  this  volume 
is  computed,  as  denoted  by  the  shaded  regions  in  the  figure 
below. 


The  initial  and  final  extensions  of  the  plug  edge  are 
special  cases.  The  initial  extension  usually  crosses  one  or 
more  cells  that  contain  projectile  as  well  as  target  material. 
In  this  case,  any  curvature  of  the  top  of  the  projectile  is 
ignored  and  the  entire  top  surface  of  the  projectile  is 
assumed  to  have  the  same  y-coordinate  as  the  beginning  point 
of  the  plug  surface.  The  partial  volume  associated  with  the 
plug,  therefore,  is  measured  from  the  approximate  position  of 
the  top  of  the  projectile  rather  than  from  the  bottom  cell 
boundary,  as  illustrated  below. 


6-15 


TAH 

Clii 

I'l.l  C. 

I'lu: 

■ MAT  I 

The  partial  volume  associated  with  the  target  is  defined  as 
the  total  cell  volume  minus  the  plug  and  projectile  partial 
volumes,  where  the  projectile  partial  volume  is  based  on  its 
mass  and  density. 

In  the  second  special  case,  when  the  plug  edge  is  ex- 
tended to  the  free  surface,  the  new  top  boundary  of  the  plug 
is  the  top  of  the  target  rather  than  a grid  line.  Ignoring 
the  small  curvature  of  the  target,  the  partial  volume  associ- 
ated with  the  plug  in  the  cell  containing  the  endpoint  is 
based  on  the  y-coordinate  of  the  new  endpoint  of  the  slipline, 
as  illustrated  below. 
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Once  the  volumes  associated  with  the  plug  material  are 
computed,  the  target  material  in  the  region  of  the  plug  ex- 
tension is  converted  into  plug  material.  The  pure  cells  to 
the  left  of  the  plug  surface  endpoint  are  converted  simply 
by  changing  the  value  of  MFLAG  for  those  cells  from  3 to  2. 
(Package  3 is  the  target;  package  2 is  the  plug.) 

The  XMASS,  SIE,  and  RHO  arrays  are  redefined  for  the 
interface  cells  that  contain  both  plug  and  target  material. 
Given  m = MFLAG(k)  - 100,  the  plug  mass  is  stored  in 
XMASS  (2,m)  and  is  defined  as  the  product  of  the  volume  of 
the  plug  material  for  that  cell  and  the  average  material 
density  for  that  cell. 


XMASS ( 2, m)  = V 


plug 


av 


(6.4) 


where 


av 


"‘cell^'^cell 
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;'he  target  material  is  stored  in  XMASS(3,m)  and  is  the 
difference  of  the  total  cell  mass  and  the  plug  mass: 

XMASS(3,m)  = AMX(k)  - XMASS(2,m)  .(6.6) 

The  density  of  both  the  plug  and  target  materials  is  defined 
to  be  p , i.e.,  RHO(2,m)  = RH0(3,m)  = p . The  specific 
internal  energy  of  both  materials  is  defined  to  be  the 
average  specific  internal  energy  for  the  cell:  SIE(2,m)  » 

SIE(3,m)  = AlX(k). 

The  procedure  for  an  interface  cell  that  will  contain 
all  three  material  packages  varies  from  that  described  above. 
In  this  case,  the  volume  of  the  target  material  is  used  to 
compute  its  mass,  so  that 


m 


liili 


XMASS(3,m)  - • p,„g„  C6.7) 


where  V^  ^=Vn-V,  -V  . *^i^.  The  density  of 

target  cell  plug  projectile  * 

the  target  material  is  used  to  define  the  plug  mass: 


XMASS(2,m)  - 


(6.8) 


This  follows  from  the  fact  that  the  plug  is  simply  a part  of 
the  target.  Therefore  the  density  and  specific  internal 
energy  of  the  plug  material  is  defined  to  be  the  same  as  the 
density  and  specific  internal  energy  of  the  target  material: 

RH0(2,m)  = RH0(3,m) 

SIF.(2,m)  = SIE(3,m) 

All  target  material  is  converted  to  plug  material  in  the 
interface  cells  to  the  left  of  the  plug  edge  and  at  the  free 
surface,  i.e.,  the  plug  quantities  are  defined: 

XMASS(2,m)  = XMASS(3,m) 

RH0(2,m)  = RH0(3,m) 

SIE(2,m)  = SIE(3,m) 

and  then  the  target  quantities  are  set  to  zero: 


XMASS(3,m)  = 0 
RH0(3,m)  « 0 
SIE(3,m)  » 0 
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DESCRIPTION  OE  INPUT  VARIABLES  AND  RESTART  PROCEDURES 

GENERAL  CAPABILITIES  AND  LIMITATIONS 

A problem  generator  which  uses  the  position  of  material 
interfaces  to  define  cell  quantities  has  been  incorporated 
into  the  HELP  code.  Currently  the  generator  includes  a sub- 
routine, TSETUP,  whicli  allows  the  user  to  easily  generate 
material  interfaces  which  are  composed  of  straight  line 
segments  or  portions  of  circles  or  ellipses.  However,  the 
method  is  applicable  to  any  topological  configuration, 
provided  the  user  generates  the  tracer  particles  wliich 
define  the  material  interfaces.  The  following  general 
capabilities  and  limitations  of  the  HELP  code  should  be 
noted  before  applying  it  to  a specific  problem: 

1.  The  .number  of  material  packages  in  a calcu- 
lation is  limited  only  by  the  dimensions  of 
the  material  arrays,  which  can  be  increased  up 
to  the  limit  of  the  user's  computer  core 
storage.  (See  Appendix  A for  dimensioning 
the  arrays.) 

2.  A single  cell  can  contain  any  number  of 
material  interfaces  and  any  number  of  distinct 
materials . 

3.  Two  material  packages  can  actually  contain  the 
same  material  (use  the  same  equation  of  state) 
but  have  different  initial  conditions 
(velocities,  energies  and/or  densities). 

4.  A package  can  be  divided  into  any  number  of 
spatially  disconnected  subpackages.  (See 
instructions  for  defining  material  tracer 
particles  in  Section  7.2.5.) 
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The  material  for  each  package  must  be 
homogeneous;  i.e.,  each  pure  cell  within 
a package  must  have  the  same  density, 
velocity  and  internal  energy.  This 
restriction  can  be  removed  by  appropriate 
user-supplied  statements  in  FILGRD. 

The  grid  boundaries  cannot  serve  as  free 
surfaces  because  the  code  assumes,  for 
example,  tliat  a package  which  extends  to 
the  right  edge  of  the  grid  is  infinite  in 
the  x-direction. 

The  grid  must  have  at  least  three  rows 
and/or  three  columns. 

The  DATA  statements  in  COMDIM  and  EQST  can 
be  redefined  if  the  user  wants  to  model  a 
material  which  is  not  in  the  EQST  list  of 
materials  modeled  by  HELP.  An  inert  material 
must  be  assigned  a material  number  from  1 
to  19,  inclusively.  A gas  or  high  explosive 
must  be  assigned  a material  number  from  20 
to  30,  inclusively.  In  order  for  the 
pressure  iteration  to  converge,  any  equation 
of  state  which  is  used  in  place  of  the 
Tillotson,  the  y-law  or  the  JWL,  must  be 
continuous  and  liave  a negative  slope  every- 
where in  the  P-V  plane,  or  at  least  prevent 
the  iteration  from  entering  any  region  which 
violates  these  restrictions.  (See  Section 
4. 6. 1.3.) 

A calculation  can  be  made  in  plane  or  in 
cylindrical  coordinates.  (See  definition  of 
IGM  in  Section  7.2.1.) 

The  top  and  right  grid  boundaries  are  trans- 

mittive.  The  left  boundary  is  reflective 

(an  axis  of  symmetry  in  cylindrical  coordinates) 
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The  bottom  boundary  is  the  only  one  which 
can  be  either  transmittive  or  reflective 
(see  definition  of  CVIS  in  Section  7.2.1). 

7.2  INSTRUCTIONS  FOR  GENERATING  A PROBLEM 

Instructions  for  generating  a liRLP  problem  are  given  in 
the  sections  that  follow.  An  abbreviated  set  of  input  in- 
structions on  keypunch  forms  is  given  in  Appendix  B. 


7.2.1  Z-Block  Variables 

Most  of  the  options,  flags,  and  problem  parameters  (such 
as  the  number  of  rows  and  columns  in  the  grid)  are  variables 
located  within  the  first  150  words  of  blank  common  and  are 
defined  by  the  first  set  of  input  cards.  These  variables 
are  called  Z-block  variables  because  they  all  follow  Z(l), 
the  first  word  of  blank  common,  and  can  be  referenced  using 
an  appropriate  subscript,  between  1 and  150,  for  Z.  For 
example,  CVIS  and  Z(27)  have  the  same  location  in  blank 
common.  Not  all  variables  in  the  Z-block  are  defined  by 
input  cards;  some  are  defined  and  updated  by  the  code  as 
the  calculation  proceeds.  All  Z-block  variables  are 
written  on  the  restart  file,  and  most  are  parameters  which 
must  be  saved  in  order  to  restart  a calculation.  However, 
the  variables  such  as  UN93  or  NUN32,  which  correspond  to 
Z(93)  and  Z(32),  respectively,  are  not  currently  being  used 
by  the  code  and  are  available  to  the  user  to  define 
additional  problem  parameters  as  the  need  arises. 

The  Z-block  variables  are  defined  by  cards  that  are 
read  by  subroutine  CARDS.  The  format  (II,  15,  II,  7E9.4)  is 
used  to  read  the  variables  lEND,  LOC,  NUMVfPC, 

(CARD(I),  1=1,  NUMWPC)  which  are  defined  as  below: 

lEND  IEND=0  indicates  that  the  variable (s)  defined 
by  this  card  are  all  typed  real . 


I 

I 


inND*2  indicates  that  the  variable (s)  defined 
by  this  card  are  all  typed  integer. 
inND“l  indicates  that  the  variable(s)  defined 
by  this  card  are  all  typed  real,  and 
this  is  the  last  card  CARDS  will  read 
until  it  is  called  again  by  INPUT. 

(Note:  a setup  deck  has  three  cards 
where  IEND=1,  a restart  deck  has  only 
two  such  cards.  Wlien  generating  a 
problem,  INPUT  calls  CARDS  three  times; 
when  restarting  a problem,  INPUT  calls 
CARDS  twice.) 

LOC  LOG  is  the  location  in  blank  common  of  the  first 
variable  being  defined  by  this  card. 

NUMWPC  NUMWPC  is  the  number  of  consecutively  stored 

Z-block  variables  of  the  same  type  being  defined 
by  this  card. 

CARD(I)  The  value  assigned  to  the  I^^’  variable.  NOTE: 

even  if  the  variable  is  typed  integer,  the  value 
must  be  entered  witli  a decimal  point.  (See 
sample  input  decks  in  Chapter  XI.) 

The  Z-block  variables  defined  by  input  cards  are  listed 
and  defined  on  the  following  pages.  The  type  (integer  or 
real)  and  location  in  blank  common  are  indicated  for  each 
variable.  The  code  assigns  non-zero  default  values  to  some 
of  these  variables  at  the  beginning  of  subroutine  INPUT. 

The  non- zero  default  values  are  given  in  parentlieses  after 
the  variable  definition.  If  no  value  is  indicated,  the 
default  value  is  assumed  to  be  zero. 


INPUT  Z- BLOCK  VARIABLES 


Column  Location 
Variable  One  in  Blank 
Name  Flag  Common 


Definition 

(Default  Values  are  in  Parentheses) 


PK(1) 


PROB 


NFRELP 


NDUMP7 


ICSTOP 


NUMREZ 


KUNITR 


Problem  number.  Any  number  from 
.0001  to  99.9999.  It  must  be 
identical  to  PROB. 

Problem  number.  (Same  value  as  PK(1).) 
Must  be  included  in  the  setup  deck. 

Gives  frequency  of  "long"  EDIT  prints, 
which  give  velocities,  energy,  com- 
pression (or  density)  and  stresses 
for  all  cells  in  active  grid.  A 
"short"  EDIT  gives  this  information 
only  for  the  first  column  of  cells, 
e.g.,  when  NFRELP®2,  every  second 
EDIT  print  is  "long,"  when  NFRELP^S, 
every  fifth  EDIT  print  is  "long",  etc. 

Gives  frequency  of  restart  tape  dumps 
relative  to  EDIT  prints.  (Program 
dumps  only  when  EDIT  prints.)  E.g., 
if  NDUMP7=5,  a restart  dump  will  be 
made  every  5th  time  EDIT  prints. 

Gives  cycle  at  which  execution  will 
stop  if  the  user  wants  to  specify 
a cycle,  rather  than  a value  of  T 
(time),  on  which  to  stop.  If  stopping 
on  time  omit  this  card. 

The  total  number  of  times  the  grid 
will  be  automatically  rezoned.  If 
not  automatically  rezoning  the  grid, 
omit  this  card. 

The  name  of  the  file  INPUT  will  read 
from.  (7) 


■W 

■ 


INPUT  Z-BL(  CK  VARIABLES 


Variable 

Name 

Column 

One 

Flag 

Location 
in  Blank 
Common 

Definition 

(Default  Values  are  in  Parentheses) 

IPR 

2 

15 

Maximum  number  of  iterations  CDT  will 
perform  when  attempting  to  equilibrate 
the  pressures  of  materials  in  a multi- 
material cell.  If  the  pressures  are 
not  within  an  epsilon  (PRCNT)  of  the 
average  pressure  after  IPR  iterations, 
an  error  exit  occurs.  (35) 

PRCNT 

16 

Convergence  criterion  tor  the  pressure 
equilibration  of  material  in  a multi- 
material cell.  All  material  pressures 
P^,  must  satisfy  the  following: 

1 (P  - P.)/P  1<PRCNT.  (.001) 

1 '■  av  1 av  ' 

KUNITW 

2 

17 

The  name  of  the  file  EDIT  and  SETUP 
will  write  on.  (7) 

I CM 

2 

21 

When  IGM  = 0 code  uses  cylindrical 
coordinates.  When  IGM  = 1,  code  uses 
plane  coordinates. 

DM  IN 

24 

The  maximum  relative  error  in  the 

energy  sum  that  the  user  wants  to 
tolerate.  The  error  is  computed  as 
follows : 


,KMAX 

\ - ETH^/ETH 

where  E.  is  total  energy  of  cell  k 
and  ETH^is  theoretical  energy  in  the 
grid- -computed  in  SETUP  and  updated  in 
HPHASE,  TPHASE,  SPHASE,  ADDENG,  REZONE, 
and  RNDOFF.  (lO'^) 

CVIS  27  A flag  that  describes  the  condition  of 

the  bottom  grid  boundary.  If  CVIS  = 0, 
the  bottom  boundary  will  be  reflective. 
If  CVIS  = - 1,  the  bottom  boundary  will 
be  transmitt ive . 
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INPUT  Z-BLOCK  VARIABLES 


Column  Location 
Variable  One  in  Blank 
Name  Flag  Common 


I MAX 


JMAX 


MAPS 


NUMSCA 


PRLIM 


Definition 

(Default  Values  are  in  Parentheses) 


The  number  of  columns  in  the  grid.  If 
planning  to  rezone  the  grid  in  the 
x-direction,  IMAX  must  be  an  even 
number.  IMAX  also  must  be  at  least  3. 

The  number  of  rows  in  the  grid.  If 
planning  to  rezone  in  the  y-direction, 
JMAX  must  be  an  even  number.  JMAX  also 
must  be  at  least  3. 

When  MAPS  > 0,  part  of  the  EDIT  print 
is  a set  of  symbolic  maps  of  the 
compression  (or  density),  pressure, 
radial  and  axial  velocities,  and 
specific  internal  energy  of  the  cells 
in  the  active  grid. 

MAPS=1  gives  a compression  map. 

MAPS=2  gives  a density  map. 

The  number  of  times  the  cycle  or  time 
interval  between  EDIT  prints  is  in- 
creased. (See  PRLIM) 

The  time  or  cycle  at  which  the  EDIT 
print  frequency  is  ciianged. -Example : 
you  wish  to  print  every  10  sec  until 

T = 10"'  sec;  and  every  10"'  sec  until 

T = 10'^  sec,  and  every  lO"^  sec  there- 
after. Set:  - 

PRDELT  =10'® 

IPCYCL  = 0. 

PRLIM  =10"^ 

PRFACT  =10. 

NUMSCA  = 2. 
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input  Z'BLOCK  variables 


Column  Location 
Variable  One  in  Blank 
'"'Same  Flag  Common 


PRDELT 


PRFACT 


II 


45 


46 


47 


12 


IPCYCL 


TSTOP 


48 


49 


50 


IPLGRT 


56 


CPe 


Definition 
fault  Values  are  in 


Parentheses) 


rncrcfsfa.b^  ?.;e 

same  factor  interval. 

omit’*NUMSCA,  PRLIM,  PRIACT. 

The  time  Csec)  '’®^^p®^j^tervals , omit 
If  printing  on  cycle  interva  , 

this  card. 

The  factor  "j!j;?a‘j%hrprint‘limit 
PlM)  rra’lSL^d.  CSee  PRUM) 

The  number  of 

have  "°"-fi”J^l2''define  the  "active" 

!S-;ifrei?r^^^ire%re%rtfvrrria., 
See-rn-l/ro  re^ocfti'erof  itr 

plus  2. 

U nf  rvcles  between  EDIT 

™!„"ts  U printing  on  time  intervals. 

omit  this  card. 

1 nf  T ftime)  at  which  the  ^ 

The  value  of  When  stopping 

calculation  will  Pj.j^g,^,QP  ^ 0),^set 

note:  This^cajd  has^a  ,1_ 

SudeHS  bfth  the  setup  de  J and  the 
restart  deck. 

The. rightmost  column  of  the  plugging 

To^i^e  ;!giror;fe^  XP  ctea 

-ara‘r’'oSit'^hifcirrwre^^rot 


4*  V>  A 


II  A.  V — — 

rvliiooinp  option. 
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INPUT  Z-BLOCK  VARIABLES 


Column  Location 
Variable  One  in  Blank 
Name  Flag  Common 


Definition 

(Default  Values  are  in  Parentheses) 


PLWMIN 


IPLGBT 


IPLGTP 


GAMMA 


NMAT 


CYCMX 


CYCPH3 


NTRACR 


The  required  minimum  specific  plastic  work  J 
(ergs/gm)  whicli  the  material  near  the  i 

top  of  the  vertical  edge  of  the  plug  | 

must  have  before  the  plug  edge  is  ex-  I 

tended  to  the  next  row.  See  Section  a 

6,2.1.  Omit  this  card  when  not  using  | 
the  plugging  option.  | 

The  bottom-most  row  of  the  plugging 
region  of  the  target.  (Usually  the  \ 

first  row  of  the  target.)  Omit  this 
card  when  not  using  the  plugging  option. 

The  top-most  row  of  the  plugging  region  j 
of  the  target.  (Usually  the  top  row 
of  the  target.)  Omit  this  card  when 
not  using  the  plugging  option.  . 

Y in  P=(y-1)pE  for  material  #20,  an  ■'< 

ideal  gas,  ^ 

i 

The  number  of  material  packages,  ex- 
cluding the  void  package. 

The  number  of  passes  through  subroutine 
INFACE  each  cycle  to  minimize  transport 
noise  near  interfaces.  (2) 

The  number  of  passes  through  subroutine 
SPHASE  (strength  phase)  each  cycle. 

When  the  calculation  is  purely  hydrodynamij 
(no  strength  effects),  set  CYCPH3*-!. (1) 

The  number  of  material  tracers  per  cell 
diagonal.  Used  by  ADDTCR  when  adding 
tracers.  ADDTCR  adds  tracers  only  when  ’ 
the  distance  between  tWo  consecutive  j 

tracers  in  a specified  region  is  .1 

greater  than  a cell  diagonal  divided  'M 
by  NTRACR.  MINX,  MAXX,  MINY,  MAXY  J 

specify  the  region  of  the  grid  considered^ 
by  ADDTCR.  NADD  specifies  the  ^ 

frequency  with  which  ADDTCR  is  called. 
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Variable 

Name 

NMXCLS 

NTPMX 

NTCC 

SIEMIN 

EMIN 

PMIN 


INPUT  Z-BLOCK  VARIABLES 


Column  Location 
One  in  Blank 
Flag  Common 


Definition 

(Default  Values  are  in  Parentheses) 


73 


78 


81 


82 


85 


86 


The  maximum  number  of  interface  cells 
to  exist  in  the  grid  on  any  one  cycle 
during  the  calculation.  This  maximum 
should  correspond  to  the  dimensions 
of  the  material  arrays  (XMASS,  SIE, 
etc.).  See  Appendix  A. 

The  maximum  number  of  tracer  particles 
per  material  package.  This  maximum 
should  correspond  to  the  dimensions 
of  the  TX  and  TY  arrays. 

When  NTCC  0,  NTCC  passive  tracer 
particles  are  initially  placed  in  the 
center  of  every  fourth  nonempty  cell 
and  subsequently  moved  through  the 
grid  with  the  material.  These  are 
used  only  for  producing  particle  plots 
and  do  not  affect  the  material  behavior. 
When  this  card  is  omitted  (NTCC  = 0)  , 
cell  centered  tracers  are  not  computed 
unless  the  plugging  option  is  being 
used,  in  which  case  they  are  automati* 
cally  generated  in  the  plugging  region 
of  the  target.  (See  Section  6.1.3) 

A cutoff  on  the  total  specific  internal 
energy  increment  of  a cell  in  TPHASE. 
This  cutoff  prevents  small  numerical 
signals  from  enlarging  the  active  grid. 
(105) 

The  minimum  value  of  specific  internal 
energy  to  be  used  in  the  ideal  gas 
equation  of  state.  (10^) 

A pressure  cutoff.  If  lP(k) | < PMIN, 
then  P(k)  =0.  (5  x 10°) 
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INPUT  Z- BLOCK  VARIABLES 


I 


Column  Location 
Variable  One  in  Blank 
Name  Flag  Common 

INTER  2 87 


NODIIMP 


NLINER 


2 105 


NVRTEX 


2 109 


Definition 

(Default  Values  are  in  Parentheses) 

A special  editing  flag  for  debugging: 

»2  for  debug  prints  in  CUT 

“3  for  debug  prints  by  EDIT  after 
HPHASE  and  SPHASE 

*5  for  debug  prints  ia  TPHASE 

*7  for  debug  prints  in  SPHASE 

=11  for  debug  prints  in  UVMOD 

■13  for  debug  prints  in  UVCALC 

(See  Section  9.2  if  more  than  one  debug 
is  desired.) 


The  flag  which  initiates  the  rezoning 
of  the  grid.  If  NUMREZ  > 0,  this  flag 
will  be  set  automatically  by  the  code. 
(See  Section  8.4)  However,  the  user  can 
force  a rezone  on  any  restart  cycle  by 
setting  REZ  ■ 1 in  the  restart  deck. 
(lEXTX  and/or  JEXTY  must  also  be  defined 
in  order  for  the  grid  to  be  rezoned.) 

When  NODUMP  = 1,  EDIT  will  not  write 
any  restart  dumps.  (This  flag  over- 
rides the  NDUMP7  option,  but  does  not 
prevent  SETUP  from  writing  the  cycle 
0 dump . ) 

The  package  number  of  the  shaped  charge 
liner  material  that  forms  the  jet.  Used 
only  when  calculating  the  collapse  of  a 
liner. 

The  second  index  of  the  void  tracer 
particle  that  is  at  the  vertex  of  the 
void  closing  region,  i.e.,  (TX(NVOID, 
NVRTEX),  TY(NV0ID,  NVRTEX))  are  the 
coordinates  of  the  vertex  point.  Used 
only  when  the  code's  automatic  void 
closing  routine  is  to  be  activated. 

(See  Section  8.4) 


' 4.V 
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INPUT  Z-BLOCK  VARIABLES 


Column  Location 
Variable  One  in  Blank 
Name  Flag  Common 


Definition 

(Default  Values  are  in  Parentheses) 


ROEPS 

110 

A round-off  epsilon  used  to  defipe 
cutoffs  in  the  calculation.  (10'“’) 

PLGOPT 

111 

If  the  user  is  generating  a plugging 
calculation,  PLGOPT  must  be  set  equal 
to  1.  The  special  plugging  mechanisms 
are  activated  only  if  PLGOPT  = 1.  (See 
Chapter  VI  before  using  the  plugging 
option. ) 

NSLD 

2 

112 

The  maximum  number  of  cells  the  user 
expects  will  contain  the  slipline  on 
any  one  cycle.  This  maximum  should 
also  depend  on  the  dimensions  of  the 
slipline  arrays  (MSLD,.  etc. ) . 

FINAL 

113 

The  final  value  of  the  stability 
fraction  used  in  determining  the  time 
step.  See  STAB.  (.4) 

NOSLIP 

2 

115 

If  no  sliplines  are  to  be  generated, 
set  NOSLIP  = 1,  which  causes  t.Le  in- 
structions and  routines  that  affect 
only  slipline  cells  to  be  skipped. 

LVISC 

2 

116 

A linear  artificial  viscosity  term 
is  added  to  cell  boundary  pressures 
if  LVISC  = 1.  (See  Section  2.2.2..  5.) 

NADD 

2 

118 

A flag  for  automatically  adding  material 
tracer  particles  in  the  region  specified 
by  MINX,  MAXX,  MIMY,  MAXY.  If  NADD  = 10 
ADD'iCR  is  called  every  cycle  which  is  a 
multiple  of  10,  (e.g.,  cycle  20,  30,  ... 
250,  260,  etc.).  NOTE:  NTRACR  (=Z(72)) 
must  also  be  defined  before  material 
tracers  can  be  added  by  ADDTCR. 

MINX 

MAXX 

2 

2 

119  ^ 

120 

These  parameters  specify  the  left  and 
right  columns  and  the  bottom  and  top 
, rows,  respectively,  of  the  region  in 

MINY 

2 

121 

which  ADDTCR  will  add  tracer  particles. 
These  variables  must  be  defined  when 

MAXY 

2 

122  , 

NADD  f 0. 

7-12 


INPUT  Z- BLOCK  VARIABLES 


Column 

Variable  One 
Name  Flag 


Location 
in  Blank 
Common 


Definition 

(Default  Values  are  in  Parentheses) 


lEXTX 


JEXTY 


STAB 


DTMIN 


CRATIO 


A rezone  flag.  The  grid  will  be  rezoned 
in  the  x-direction  only  when  lEXTX  « 1. 
Omit  this  card  if  the  grid  is  not  to  be 
rezoned,  or  if  it  is  to  be  rezoned  in 
the  y-direction  only. 

A rezone  flag.  The  grid  will  be  rezoned 
in  the  y-direction  only  when  JEXTY  * 1. 
Omit  this  card  if  the  grid  is  not  to  be 
rezoned,  or  if  it  is  to  be  rezoned  in 
the  x-direction  only. 

The  initial  value  of  the  stability 
fraction  used  in  determining  the  time 
step.  If  FINAL  * 0,  STAB  is  constant. 
Otherwise,  its  value  is  doubled  each 
cycle  until  it  reaches  the  value  of 
FINAL.  (10’3) 

The  minimum  value  for  the  time  step. 

The  program  gives  an  error  exit  if  CDT 
calculates  a At  < DTMIN.  (10"^^) 

The  compression  of  the  material  in 
adjacent  cells  is  computed.  If  the 
ratio  of  these  cells’  compressions  is 
greater  than  CRATIO,  their  pre'.sures 
are  inverse  compression  weighted 
and  their  velocities  are  compression 
weighted  to  define  cell  boundary  values 
in  HPMASE.  (See  Section  2 . 2 . 2 . 4 . ) (10^) 


BBAR 


EMOB 


A constant  used  in  the  calculation  of 
the  local  sound  speed  of  materials 
other  than  ideal  gases  and  high  explo- 
sives. C = C^  BBAR  • >/|P|  . (.5) 


Dummy  end  card  of  a setup  deck.  This 
card  is  read  by  CARD^  after  all  other 
input  cards  have  been  read.  “TlTis  card 
must  not  be  included  in  a restart  deck. 
Set  EMoB  = 0 in  the  setup  deck. 
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7.2.2  Cell  Dimensions 


The  second  set  of  inpjit  cards  define  the  cell  dimensions. 
The  X and  y dimensions  of  the  cells  can  vary;  however,  the 
dimensions  of  contiguous  cells  should  not  vary  by  more  than 
10%  except  in  parts  of  tlie  grid  well  beyond  the  region  of 
interest.  Large  aspect  ratios  (e.g.,  AX/AY  > 2)  are  also 
discouraged  in  significant  regions  of  the  grid. 

The  cards  which  define  the  DX,  DY  arrays  are  road  by 
subroutine  SETUP.  A set  of  cards  defining  DX  are  read,  then 
a set  defining  DY.  Each  card  specifies  up  to  four  different 
DX  or  DY  values.  The  informatioa  is  read  into  the  NNT  and 
TEMP  arrays:  (NNT(L),L  = 1,4),  (TEMP(L),L  = 1,4),  where  NNT(L) 
is  the  number  of  columns  (or  rows)  that  have  a DX  (or  DY)  dimen- 
sion equal  to  TEMP(L).  After  the  x-dimension  (DX)  of  all 
columns  has  been  defined,  the  next  NNT  is  set  equal  to  999, 
which  indicates  that  the  entire  DX  array  is  defined.  Then 
the  DY  array  is  similarly  defined  by  the  next  group  of  cards. 
After  the  y-dimension  (DY)  of  all  rows  has  been  defined,  the 
next  NNT  is  set  equal  to  999,  indicating  that  all  y-dimensions 
have  been  defined.  An  error  exit  will  result  if  the  user  does 
not  define  exactly  IMAX  DX's  and  exactly  JMAX  DY's. 

7.2.3  Initial  Density,  Velocity,  and  Specific  Internal 
Hnergy  of  Each  Material  Package 

In  addition  to  the  material  tracer  particle  positions,  the 
definition  of  material  packages  requires  identification  of 
MFLAG(2),  the  flag  of  the  cell  in  the  bottom  left  corner  of 
the  grid,  and  the  material  code  number  and  initial  conditions 
of  each  material  package. 

The  card  which  follows  those  defining  the  cell  dimensions 
defines  MFLAG(2).  The  user  should  be  able  to  predict  from  the 
placement  of  the  interfaces  in  the  grid  whether  cell  k = 2 
will  be  a pure  cell,  an  interface  cell,  or  a void  cell  at 
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t = 0.  If  it  will  be  a pure  cell  of  package  n,  then  set 
MFLAG(2)  =»  n.  If  it  will  be  an  interface  cell,  then  set 
MFLAG(2)  ■ -1.  (In  the  latter  case  its  value  will  be 
changed  to  be  greater  than  100  after  the  interfaces  are 
processed  by  CALFRC  and  VOLFND) . If  it  will  be  a void 
cell,  then  set  MFLAG(2)  » 0. 

After  MFLAG(2)  is  defined,  a set  of  cards  which  defines 
the  material  code  number  and  initial  conditions  of  each 
material  package  is  read.  The  first  card  in  this  set  defines 
material  package  1,  the  second  card  defines  material  package 
2,  etc.  NMAT  cards  must  therefore  be  read  to  define  all  the 
material  packages. 

The  material  code  numbers  for  24  materials  are  listed 
in  EQST  as  well  as  in  Tables  2.1  and  2.2  of  this  report. 

The  material  code  numbers  are  stored  in  the  MAT  array  and 
are  used  whenever  the  equation  of  state  constants  (Pq»A,B,y» 
etc.)  are  referenced.  These  equation  of  state  constants  are 
defined  in  DATA  statements  in  EQST  and  in  COMDIM,  the 
"included"  element.  For  example,  if  package  3 is  lead 
(material  code  number  10)  then  MAT(3)  = 10,  ESCAPA(IO)  is 
its  bulk  modulus.  A,  and  RHOZ(IO)  is  its  normal  density,  p^. 

The  density  (g/cc),  specific  internal  energy  (ergs/g) 
and  the  radial  and  axial  velocities  (cm/sec)  of  all  the 
material  in  each  package  are  specified  by  these  cards.  The 
knowledgeable  user  can  vary  these  quantities  within  the 
package  by  adding  to  subroutine  FILGRD  coding  which  would 
redefine  the  appropriate  variables. 

When  defining  the  initial  density,  RHOIN,  of  a material, 
the  user  should  be  aware  of  the  normal  density,  RHOZ, 
specified  by  the  code  fcr  that  material.  (See  Table  2.1  and 
2.2.)  The  user  can  change  the  definition  of  for  a given 
material  by  redefining  the  appropriate  RHOZ  variable  in  the 
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DATA  statement  in  COMDIM, the  "included"  element.  The  user 
can  also  specify  an  initial  density,  RHOIN,  which  is  dif- 
ferent from  a material's  normal  density,  RHOZ. 

See  Appendix  B for  the  exact  format  of  the  cards  which 
define  MI-’LAG(2)  and  tlie  material  code  numbers  and  initial 
conditions  of  the  material  packages. 

7.2.4  Strength  Constants  for  Each  Material  Package 

The  shear  yield  strength  and  the  tensile  strength  of 

each  material  package  are  defined  by  the  next  set  of  cards. 

The  shear  yield  strength  constants,  Y , Y, , Y^,  E , G,  are 

o o l ^ m 

all  set  to  zero  if  the  material  has  no  yield  strength.  If, 
however,  the  material  has  no  tensile  strength,  AMDM(n)  is 
set  to  1.  The  definition  of  the  yield  strength  terms  is 
given  in  Section  2. 3. 2. 3.  The  tensile  failure  threshold  is 
described  in  Section  2.3.3. 

The  setup  deck  must  contain  NMAT  cards  defining  strength 
properties,  even  when  none  of  the  packages  have  strength  and 
SPHASE  is  being  hypasse<!. 

See  Appendix  B for  the  exact  format  of  the  cards  which 
define  the  strength  constants. 

7.2.5  Material  Tracer  Particles 

A material  package  boundary  is  generated  by  dividing  it 
into  a series  of  straight  lines,  and/or  arcs  of  circles  and/or 
arcs  of  ellipses. 

A set  of  2 or  3 cards  is  read  by  subroutine  TSETUP  to 
define  each  one  of  these  segments.  The  first  of  these  cards 
specifies:  (1)  LTYPE,  what  kind  of  a segment  to  generate  - a 

horizontal  line,  a vertical  line,  a diagonal  line,  an  arc  of 
a circle,  or  an  arc  of  an  ellipse;  (2)  MPN,  which  material 
package  is  being  generated;  and  (3)  NPTS,  how  many  tracer 
particles  to  put  along  the  segment.  If  the  segment  is  a 
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straigljt  line,  only  one  other  card,  which  defines  the  end- 
points of  the  line  in  centimeters,  is  read.  However,  if 
the  segment  is  an  arc,  two  other  cards  are  read  which  define 
the  beginning  and  ending  angles  of  the  arc  in  radians  and 
the  location  of  its  center,  or  foci,  in  centimeters. 

iVhen  defining  the  last  segment  of  a material  package  or 

* 

subpackage  boundary,  the  user  must  assign  a negative  value 
to  LTYPE.  This  signals  TSETUP  to  generate  a dummy  tracer 
particle  which  marks  the  end  of  a tracer  string  for  a sub- 
package or  package.  The  coordinates  of  this  dummy  particle 
are  (-1000,  0),  and  an  error  exit  will  result  if  this  dummy 
particle  is  not  the  last  particle  in  the  set  of  tracers  for 
a material  package.  After  all  segments  of  all  packages,  in- 
cluding the  void  package,  have  been  generated,  a flag  card, 
setting  LTYPE  * 100,  is  read  by  TSP-TUP.  This  signals  that 
all  the  material  tracer  particles  have  been  defined. 

As  the  segments  are  generated  they  must  form  a continuous 
boundary  interrupted  only  by  the  grid  boundaries.  (It  is 
unnecessary  to  generate  tracers  along  the  grid  boundaries.) 
Furthermore,  the  segments  need  to  be  generated  in  an  order 
such  that  the  interior  of  the  package  lies  on  the  left  as 
one  moves  between  any  consecutive  pair  of  tracers.  This 
latter  requirement  also  dictates  which  endpoint  must  be 
specified  first  on  the  input  cards. 

Since  every  segment  of  an  interface  is  described  by 
two  identical  strings  of  tracer  particles  (one  for  each 
material  package  lying  on  opposite  sides  of  the  interface) , 


A material  package  can  be  divided  into  any  number  of 
spatially  disconnected  subpackages.  The  tracer  string 
for  each  subpackage  must  end  with  a dummy  tracer. 
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the  number  of  tracers  generated  along  a segment  must  be 
identical  for  both  packages,  the  order  of  the  two  sets  of 
tracers  along  that  segment  being  exactly  opposite. 

A good  rule  to  use  in  deciding  how  many  tracers  to 
place  along  a segment  is  to  generate  two  per  cell  edge. 

In  regions  where  great  distortions  are  likely  to  occur,  the 
tracers  may  be  generated  more  densely,  although  activating 
ADDTCR  (see  definition  of  NADD  in  Z-block  variables)  helps 
to  insure  that  the  tracers  will  not  become  too  sparse. 

The  user  cannot  generate  more  than  NTPMX  (a  Z-block 
input  variable)  tracers  for  a given  material  package 
boundary;  in  that  event  TSETUP  will  print  a message  and 
stop.  The  definition  of  NTPMX  should  match  the  dimensions 
of  the  tracer  particle  arrays,  TX,  TY. 

See  Appendix  B for  the  exact  format  of  the  cards  which 
define  the  material  tracer  particles. 


7.2.6  Slipline  Endpoints 

When  sliplines  are  being  generated,  SETUP  reads  a set 
of  cards  which  defines  the  siipline  endpoints  and  specifies 
each  material  package  as  a master,  a slave,  or  neither,  (If 
there  are  no  sliplines  in  the  calculation  and  if  NOSLIP 
(a  Z-block  input  variable)  is  set  equal  to  1,  the  slipline 
cards  will  not  be  read  by  SETUP  and  must  be  omitted  from 
the  input  deck.) 

One  card  is  read  for  each  package;  the  first  card 
defines  the  first  package,  the  second  card  defines  the 
second  package,  etc.  Each  of  the  first  two  variables, 
MASTRD,  NSLAVD,  is  set  either  to  the  package  number  or  to 
zero,  depending  on  whetlier  the  package  is  a master  or  slave. 
Each  of  the  next  two  variables  NBGMD,  NBGSD  is  set  either  to 
the  index  of  the  first  tracer  of  the  slipline  or  to  zero, 


7-18 


t 


again  depending  on  whether  the  ^^NDSD,  is 

^:rri::;r^  ti:  rder:i  the  iast  tracer  ei  the  Siipune 
::  to  tore,  depending  on  whether  the  package  ts  a .aster 

slave . . 

:::  rr::.:: 

irster'^Tt'o::' nrt^rtfe":hi:h";:c'rages  are  the  .asters 
Id  which  are  the  slaves.^  cellTal  centrinriire 

:ran^:::”is::;a::::the-^^^^^^^^^^^^^ 

rac'rgL'u  T^Z'Z7Z7t^..  sa.e  velocity  co.ponents. 

'"’“see  Appendix  B for  the  exact  for.at  of  the  cards  which 
define  the  slipline  endpoints. 

7 2 7 ..Cinition  of  HigB  Pvpiosive  hetonationloints 

SPTUP  reads  two  sets  of  cards  in  order  to  define  the 
P7ia.ary  and  secondary  initiation  points  of  one  or  .ore  g^ 

:-rorrr;J-::rr,, 
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starting  with  the  point  having  the  minimum  detonation  time. 
There  is  one  card  for  each  detonation  point,  plus  an  end 
flag  card. 


The  second  set  of  cards  defines  the  approximate  region 
of  detonation  associated  with  each  detonation  point.  These 
regions  can  overlap.  There  is  one  card  in  this  set  for  eacli 
detonation  point,  hut  no  end  flag  card. 

If  the  calculation  does  not  involve  a high  explosive, 
one  blank  card  must  be  inserted  in  the  deck. 

See  Appendix  B for  the  exact  format  of  the  cards  which 
define  the  detonation  points. 


7.3  RESTART  PROCEDURES 

During  the  calculation,  HELP  periodically  writes  on  a 
tape  or  drum  file  tlie  problem  parameters  and  the  current 
state  of  the  material  in  each  cell.  By  reading  this  file 
and  a few  input  cards,  the  code  can  "restart"  and  continue 
a calculation  from  an  intermediate  point. 

In  general,  a restart  deck  consists  of  five  parts: 

1.  the  heading  card 

2.  the  "restart"  card 

3.  cards  redefining  various  Z-block  variables 
(if  desired) 

4.  the  "end  of  data"  card 

5.  cards  defining  slipline  endpoints  (if  the 
slipline  option  is  activated;  i.e.,  if 
NOSLIP  = 0) . 

The  heading  card  can  contain  any  alphanumeric  symbols  between 
columns  2 and  72.  The  "restart"  card  defines  three  variables 
in  the  PK  array  whicli  are  described  below.  (Note:  since  the 
"restart"  card  contains  three  words,  it  must  have  a "3"  in 
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column  7.  See  Section  7.2.1.)  Also,  since  it  is  the  only 
card  processed  by  CARDS  before  the  restart  file  is  read,  it 
must  have  a "1"  in  column  1.  The  third  part  of  the  restart 
deck  allows  the  user  to  redefine  any  of  the  Z-block 
variables  he  desires;  e.g.  variables  which  make  the  code 
add  tracer  particles,  change  edit  frequency,  obtain  debug 
prints,  etc.  If  the  user  is  stopping  the  calculation  on 
cycles,  he  must  define  a new  stop  cycle,  ICSTOP,  and  this 
is  the  only  case  in  which  there  must  be  a card  in  this 
third  part  of  the  restart  deck.  Tlie  fourth  part  of  the 
restart  deck  is  the  "end  of  data"  card.  This  card,  which 
must  have  a "1"  in  column  one,  redefines  TSTOP,  the  value  of 
time  at  which  the  current  run  will  stop.  If  the  calculation 
is  stopping  on  cycles,  rather  than  time,  a zero  value  is 
entered.  The  last  part  of  the  restart  deck,  which  defines 
the  slipline  endpoints,  is  necessary  only  if  the  slipline 
fcapability  is  being  used  in  the  calculation.  (See  Section 
7.2.6  for  discussion  of  defining  the  slipline  endpoints.) 

As  noted  in  Section  8.3,  the  slipline  can  be  completely 
redefined  or  even  removed  during  a restart.  The  user  should 
be  particularly  aware  of  the  possible  change  of  endpoint  in- 
dices if  ADDTCR  (see  Section  8.2)  and/or  VDCLOS  (see  Section 
8.4)  have  been  called  during  the  previous  run.  In  general, 
unless  the  user  wishes  to  redefine  the  endpoints  of  the 
slipline,  he  sl\ould  input  those  indices  listed  in  the  EDIT 
print  of  the  cycle  from  which  the  calculation  is  being 
restarted . 


The  restart  deck  of  6 cards  listed  below  could  be  used 
to  restart  at  cycle  78  a calculatipn  whose  problem  number 
is  26.5  and  to  run  it  until  T = 3.0  ysec.  A "short"  EDIT 
print  of  the  restart  cycle  is  indicated  by  PK(3)  * — 2. 

The  third  card  redefines  NADD,  which  controls  the  frequency 
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of  the  calls  to  subroutine  ADDTCR  (let's  assume  that  MINX, 
MAXX,  MINY,  MAXY  were  defined  previously).  Until  NADI)  is 
redefined  again  ADDTCR  will  be  called  every  cycle  which  is 
a multiple  of  10.  The  slipline  cards  indicate  package  1 is 
a slave  and  package  2 is  a master;  23  tracers  in  each 
package  define  the  slipline. 


j 8-16  1 17-25  I 26-34  | 35-43  | 44-52  j 53-61  j 62-70  | 

I 2 3‘)5a7»’o  I2  3*)567^<»{j(  . | 2 H 1 23jM54789o  J ^3‘)S67 B vc  i|23 MSA 789o 


tBl|324.! 
I lali  to. 

C 1 

2 n' 


IHf>n  - aLUiiIuUn  P‘PACT 

79.  «2. 


-CA 

n 

147 


Let's  assume  that  this  same  problem  was  to  be  restarted 
on  cycle  325,  the  slipline  was  to  be  removed  and  the  calcula- 
tion was  to  be  stopped  on  cycle  326.  The  following  5 input 
cards  would  be  used; 


! 8-16  I 17-25  1 26-34  j 35-43  I 44-52  1 53-61  I 62-70  1 

1 2 3 H b 4 7 8 9 7 I 2 3 9^4  7 8 9j  I 2 3 J 54  7 0 9 ^ | 2 3 9 b 47  8 9 p I 2 3 9 b>7  7 8 9 q 


1 151324,5 

2 I ) bl  I • 

2 7132a. 

I Solo. 


iKnn  - /LUnlHUi  Impact 
325.  -2. 


Card  3 sets  NOSLIP  = 1,  which  deactivates  the  slipline 
option.  Thus  it  is  not  necessary  to  define  the  slipline 
endpoints  in  this  restart  deck. 


n 
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rt  should  be  noted  that  any  blank  common  variable  that 
is  not  in  the  Z-block  can  be  changed  using  the  CARDS 
routine  format  if  the  user  can  compute  its  location  in  blank 
common.  However,  usually  it  is  easier  to  code  the  changes 
into  subroutine  INPUT  (following  statement  number  40)  after 
the  restart  tape  has  been  read  and  the  pressure  array  has 
been  initialized.  Of  course,  this  is  the  only  way  variables 
not  in  blank  common  can  be  redefined  by  the  user. 


In  the  course  of  some  calculations,  problems  can  arise 

because  of  the  position  of  the  material  tracer  particles. 

% 

This  is  most  likely  to  occur  in  calculations  which  involve 
severe  velocity  gradients  or  which  employ  VDCLOS,  the  auto- 
matic void  closing  routine  (see  Section  8.4).  In  the  first 
case,  the  speed  at  which  the  tracer  particles  become  sparse 
makes  it  easier  for  the  tracers  to  cross  each  other, 
destroying  the  logic  of  FRACS.  In  the  second  case,  it  may 
happen,  for  example,  that  a tracer  particle  may  pass  to  the 
left  of  the  vertex  point.  This  could  create  a situation  in 
which  VDCLOS  forces  the  tracer  string  to  cross  itself  when 
the  tracers  are  moved  to  the  new  vertex  position. 

If  these  situations  arise,  it  is  necessary  for  the  user 
to  add,  delete,  or  reposition  some  of  the  interface  tracer 
particles  on  a restart  cycle  (probably  the  last  one  whic*^  was 
written  before  the  problem  arose)  in  order  for  the  calculation 
to  proceed  smoothly.  The  procedure  is  not  difficult,  but.  must 
be  carefully  thought  out  before  being  implemented.  Basically 
these  are  four  principles  which  must  be  firmly  adhered  to: 

1.  Every  interface  tracer  of  each  package,  including 
the  void,  has  a corresponding  tracer  from  the  other 

t 

package(s)  which  adjoins  that  interface.  In  the 
case  of  a corner  point  among  more  than  two  packages. 
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such  as  the  intersection  of  the  liner,  casing,  and 
HE  package  in  a shaped  charge  calculation,  there  may 
be  more  than  one  other  corresponding  tracer.  If  any 
tracers  from  one  package  are  added,  deleted,  or 
moved  to  a new  location,  its  corresponding  twin(s) 
must  also  be  added,  deleted  or  moved  to  the  identical 
new  location. 

2.  If  tracers  are  being  added  to  or  deleted  from  a 
material  or  void  package  n,  the  value  of  NMP(n), 
which  stores  the  number  of  particles,  including  the 
dummy  endpoint,  defining  the  package  n boundary, 
must  be  increased  or  decreased  accordingly. 

3.  When  adding  or  deleting  tracers,  the  user  must  in- 
sure that  the  dummy  endpoint  (-1000,  0)  is  the 
value  of  the  final  tracer  of  each  package  and  sub- 
package . 

4.  When  adding  or  deleting  tracers  from  calculations 
which  use  the  automatic  void  closing  routine, 
and/or  the  slipline  option  of  the  code,  care  must 
be  taken  to  insure  that  the  values  of  NVRTEX  and/or 
the  NBGMD,  NENUMD,  NBGSD , NHNDSD  arrays  reflect  the 
changes  being  made. 

An  example  of  a calculation  which  uses  the  automatic  void 
closing  routine  will  serve  to  illustrate  the  above  procedures. 
Consider  the  situation  in  Figure  7.1,  in  which  the  package  2 
tracer  particles  near  the  vertex  of  the  void  closing  region 
have  become  too  sparse  to  satisfy  the  conditions  considered 
by  VDCLOS,  i.e.,  even  if  the  pressure  is  positive  in  the  cell 
containing  the  vertex  point,  the  angle  between  the  two 
surfaces  is  greater  than  .2  radians  and  is  becoming  larger. 
(See  Section  8.4.2.)  Therefore,  this  void  must  be  closed  by 
t)ie  user. 
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In  this  example,  the  current  vertex  of  the  void  closing 
region  is  defined  by  the  79^  package  1 particle,  the  17^ 
package  2 particle  and  the  140^^  free  surface  particle 
(package  3).  Let  us  assume  the  user  decides  the  new  vertex 

i"  K 

point  will  be  defined  by  the  18th  package  2 tracer.  The  75^ 
package  1 tracer  will  be  moved  to  the  new  vertex  point,  while 
package  1 tracer  particles  76  through  78  and  free  surface  » 

tracers  140  through  144  will  be  removed.  Figure  7.2  illus- 
trates the  modified  interface  positions  resulting  from  these 
changes  to  the  tracer  particles.  Note  that  the  139^^  free 
surface  tracer  is  now  at  the  vertex,  so  NVRTEX  must  be  set 
equal  to  139  for  VDCLOS  to  operate  correctly. 

These  changes  to  the  tracer  particles  can  be  made  in 
subroutine  INPUT  as  indicated  by  the  following  lines  of  FORTRAN: 


i 

I 

I 

I 

I 

I 


c 

30 

MO 

c 

c 

r 

C 


C 


351 


C 


352 


C 


•••  initialise  T'-STORAgE  # 

DO  MC 
P ( M «C  *C 


• ••  rHA  luL  COHtiOij 

I F ( NC  • nT  . I ‘K, » (,'J  TO  355 

TX(  I ,7S)«TX(2,UI 
Ty(i,75)»TY«2,lB> 

LO  3bl  L»7V,96 
T X < I ,L-3  I »T  X ( I tL  I 

Ty(l,L-3)«lYU»Ll 
CONTI  NUf. 

NMP  I I ) •'»3 

DO  352  L-IlSiZlV 

TXC3,L-5)«Tx(3tu> 

TY(3,l.-5»«TY(3»U 

continue 

NMP ( 3 ) *2  I M 


yapiaRLEs  here  on  a 


lines  added 
to  INPUT  to 
y change  tracer 
particle 
positions 


c 

355 


NVR  ’EX»139 
CONTINUE 


V 


ke5Tap  r 


I..  - 


7-26 


I (grid  lines) 

Figure  7. 2 --Tracer  particle  positions  after  being  redefined  by  user. 


In  this  cast,  the  changes  arc  being  made  after  cycle 
In  this  , restart 

145  has  bcfin  read  off  the  restar 

dump  will  have  the  unmodified  positions  (Figure  7.1)  of  t 

^0^  rticles;  however,  in  the  next  cycle  to  be  exocu  ed, 

C tnfaCE  will  process  the  modified  tracers  as  they 

number  x4b,  INFACE  win  pro^c  

appear  iu  Figure  7.2,  and  the  next  restart  dump  wrll 
the  modified  tracer  positions. 
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INTERACTIVE  CAPABILITIES 

There  are  several  ways  the  user  can  interact  with  the 
calculation  as  it  progresses.  The  area  encompassed  by  the 
grid  can  be  enlarged  and  the  resolution  reduced  by  rezoning; 
the  density  of  the  tracer  particles  in  a specified  region 
can  be  increased;  interfaces  can  become  slip  surfaces  or  slip 
surfaces  can  be  removed;  and  a vertex  of  two  intersecting 
free  surfaces  can  be  specified  where  the  code  will  automat- 
ically "close  a void."  These  interactive  capabilities  will 
be  discussed  in  more  detail  in  the  sections  that  follow. 

8.1  RE ZONING  THE  GRID 

The  HELP  rezone  reduces  the  resolution  of  the  grid  by 
combining  cells  and  enlarging  the  area  encompassed  by  the 
grid,  keeping  the  number  of  cells  in  the  grid  constant.  The 
HELP  rezone  combines  two  adjacent  cells  into  one.  Therefore, 
the  grid  can  be  rezoned  radially  and/or  axially  only  if 
there  is  an  even  number  of  cells  in  the  direction  to  be 
rezoned.  The  user  should  note  that  this  rezone  routine  is 
not  capable  of  rezoning  plugging  problems. 

8.1.1  Combining  Cells 

When  two  cells,  k and  kn,  are  combined  into  one,  the 
masses  of  the  two  cells  are  summed,  and  the  velocity 
components  (u',v*)  of  the  new  cell  are  a mass  weighted 
average  of  the  velocity  components  of  the  original  cells, 
thereby  conserving  mass  and  momentum. 
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m'  = m(k)  + m(kn) 

u'  = (u(k)m(k)  + u(kn)m(kn))/(m(k)  + m(kn)) 

v'  = (v(k)m(k)  + v(kn)m(kn))/(m(k)  + m(kn)) 

The  material  masses  and  velocities  of  interface  cells  are 
averaged  in  the  same  manner.  For  example,  given  m=MFLAG (k) - 100 
and  mn*MFLAG(kn) -100,  the  mass  of  material  n in  the  new  cell, 
will  be 

m^'^  = XMASS(n,m)  + XMASS(n,mn) 

The  internal  energy  of  the  new  cell  is  defined  so  as  to 
conserve  the  total  energy  of  the  two  cells  being  combined. 
Therefore,  the  internal  energy  of  the  new  cell,  IE*,  becomes 
the  sum  of  the  internal  energies  of  the  two  cells  plus  the 
thermalized  kinetic  energy  resulting  from  the  definition  of 
the  new  cell  velocity  components. 

IE*  = [lE(k)*mrkj  + IE(kn)*m(kn)  + 

m(k) • [(u(k)-u')^  + (v(k)-v’)^]  + 

m(kn)  • [ (u(kn) -u* ) + (v(kn) - v* ) ^]] / (m (k)  + m(kn)) 

The  internal  energies  of  the  materials  in  interface  cells 
are  redefined  in  the  same  manner. 

The  deviator  stresses  of  the  new  cell  are  a mass  weighted 
average  of  the  deviator  stresses  of  the  two  cells,  whereas 
the  slipline  angle  and  the  detonation  time  associated  with 
the  new  cell  are  simple  averages  of  those  quantities  in  the 
two  original  cells. 

8.1.2  Adding  Material 

The  code  automatically  adds  material  above  and/or  to 
the  right  of  the  original  grid,  depending  upon  the  direction 
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in  which  the  grid  is  being  rezoned.  If,  for  example,  the 
grid  is  rezoned  in  the  x-direction,  IMAX/2  cells  are  added 
to  each  row,  since  the  combination  of  pairs  of  original 
cells  has  resulted  in  the  original  area  being  encompassed 
by  half  as  many  cells.  The  cells  that  are  added  on  the 
right  have  the  same  x-dimension  as  the  cell  in  the  last 
column  of  the  original  grid.  These  added  cells  are  also 
assumed  to  have  the  same  initial  properties  (density, 
velocity  and  internal  energy)  that  the  material  had  in  the 
initial  problem  setup. 

If  the  last  coll  in  a row  is  an  interface  cell,  the 
code  assumes  that  the  interface  extends  throughout  the  row 
being  added;  therefore,  all  the  added  cells  ift  that  row 
will  be  interface  cells.  The  user,  in  this  case,  should 
have  generated  tracer  particles  on  a horizontal  line  beyond 
the  area  of  the  original  grid.  (Note:  The  rezone  assumes 
that  interfaces  extending  beyond  the  original  grid  are 
horizontal  in  the  x-direction  and  vertical  in  the  y-direction. 
Also,  the  rezone  will  not  automatically  generate  a free- 
surface  and  a void  unless  they  are  an  extension  of  a free 
surface  and  a void  in  the  original  grid.)  Furthermore,  the 
rezone  routine,  as  it  exists,  does  not  add  material  on  the 
left  (in  plane  coordinates)  nor  below  the  grid.  However, 
the  knowledgeable  user  can  remove  these  limitations  quite 
easily. 


8.1.3  Activating  the  Rezone 

Basically  the  code  implements  two  methods  for  rezoning 
the  grid.  One  of  these  is  an  automatic  rezone,  in  which  the 
code  automatically  determines,  during  a run,  when  to  rezone, 
does  so,  and  continues  the  run  to  the  desired  final  time  or 


8-3 


cycle.  The  other  method  is  a rezone  which  is  triggered  by 
the  user  through  the  restart  deck.  No  matter  which  method 
is  being  used,  the  user  must  set  the  Z-block  variables 
lEXTX  and/or  JEXTY  to  1 to  indicate  the  direct  ion (s) , radial 
and/or  axial,  respectively,  in  which  the  grid  is  to  be 
rezoned.  As  noted  before,  the  grid  dimension(s)  (IMAX  and/or 
JMAX)  in  the  direction(s)  to  be  rezoned  must  be  an  even 
number. 

8. 1.3.1  Automatic  Rezone 

In  order  to  activate  the  automatic  rezone  capability, 
the  user  must  define  the  Z-block  variable  NUMREZ  to  be  the 
maximum  number  of  automatic  rezones  he  wishes  to  allow. 
Subroutine  TPHASE  then  tests  to  see  if  any  mass  is  trans- 
ported through  the  right  grid  boundary  (if  lEXTX  = 1)  and/or 
the  top  grid  boundary  (if  JEXTY  =1).  If  any  mass  is  lost 
through  the  boundary  being  tested,  the  rezone  flag  REZ  is 
set  to  1 and  the  code  rezones  on  the  next  cycle. 

The  code  automatically  edits  both  before  and  after  the 
rezone.  Also  the  value  of  NUMREZ  is  decreased  by  1 follow- 
ing each  rezone  so  that  a maximum  of  NUMREZ  automatic 
rezones  is  allowed. 

8. 1.3. 2 User-activated  Rezone 

The  user  has  the  option  of  triggering  the  rezone 
through  the  restart  deck  by  defining  the  rezone  flag 
REZ  * 1.  In  this  case  the  pickup  cycle  is  rezoned  in  the 
directions  indicated  by  the  lEXTX  and  JEXTY  flags.  The 
code  gives  an  edit  following  the  rezone  and  continues  the 
execution  of  the  problem. 
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8.2  AUTOMATIC  ADDITION  OF  MATERIAL  TRACHR  PARTICLES 


By  defining  the  appropriate  variables  (NADI),  NTRACR, 
MAXX,  MINX,  MAXY,  MINY) , the  user  can  insure  that  the 
tracer  particles  will  not  become  too  sparse  in  the  regions 
of  the  grid  where  there  are  great  expansions  or  distortions 
of  the  material  interfaces. 

Every  cycle  which  is  a multiple  of  NADD,  subroutine 
ADDTCR  is  called.  ADDTCR  searches  the  region  bounded  by 
the  MINX,  MAXX  codumns  and  the  MINY,  MAXY  rows.  If  a pair 
of  consecutive  tra^rs  in  that  region  are  more  than 
1/NTRACR  of  a cell  diagonal  apart,  ADDTCR  adds  enough 
tracers  between  the  pair  to  achieve  the  desired  density  of 
NTRACR  tracers  per  cell  diagonal.  (Section  7.4  describes 
how  tracers  can  be  "manually"  added  by  inserting  coding 
changes  in  INPUT  when  a calculation  is  restarted.) 

If,  by  adding  a set  of  tracers,  tlie  dimension  (NTPMX) 
of  the  tracer  particle  arrays  would  be  exceeded,  ADDTCR 
instead  prints  a warning  and  sets  NADD  to  zero.  The  code 
will  not  call  ADDTCR  until  NADD  is  redefined  in  the  restart 
input  deck,  which,  in  this  case,  should  be  done  only  after 
the  tracer  particle  arrays  have  been  enlarged  and  the  code 
recompiled.  Also,  ADDTCR  redefines  the  slipline  endpoint 
indices  and  the  void  closing  vertex,  NVRTEX,  when  necessary. 

8 . 3 REDEFINING  SLIPLINES 

The  sliplines  in  HELP  can  be  lengthened,  shortened, 
created  or  eliminated  each  time  a calculation  is  restarted, 
simply  by  redefining  the  endpoints  of  the  slipline.  During 
the  first  cycle  following  a restart  the  code  will  reflect 
the  change  in  the  slipline  definitions  when  it  computes  the 
value  of  THETA(m)  for  each  interface  cell. 


An  interface  cell  whose  interface  is  changed  from  a slip 
to  a non-slip  surface,  will  have  all  material  velocity  component 
set  equal  to  the  cell  centered  velocity  components  and  the 
resulting  thermalized  kinetic  energy  will  be  added  to  the 
material  internal  energies.  This  will  be  done  automatically 
in  TPIIASE. 

The  material  velocity  components  of  interface  cells 
that  have  changed  from  non-slip  to  slip  cells,  will  gradually 
become  different  from  the  cell -centered  values  as  the  master 
and  slave  materials  are  treated  differently  in  TPHASE  and 
HPHASE. 

It  should  be  remembered  (when  redefining  the  endpoints 
of  the  slipline),  that  a slip  cell  has  no  shear  yield 
strength. 

8.4  CLOSING  A VOID 

The  continuous  velocity  field  used  to  move  tracer  parti- 
cles in  HELP  prevents  two  particles  on  a collision  course 
from  ever  actually  meeting.  As  the  particles  come  closer  and 
closer  to  one  another  their  velocity  fields  become  more  and 
more  alike  until  they  are  being  moved  with  almost  identical 
velocity  components.  Therefore,  if  two  free  surfaces  are 
moving  toward  one  another,  they  will  never  quite  meet,  and 
the  void  will  never  quite  close. 

The  void  closing  routine  in  HELP  will  automatically 
close  the  void  in  one  specified  region  of  the  grid  where 
two  free  surfaces  initially  have  a common  point;  this  common 
point  is  called  the  vertex  in  the  discussion  that  follows. 
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8.4.1  Identifying  Void  Closing  Region 

The  vertex  point  is  identified  by  an  input  variable, 
NVRTEX,  where  (TX(NVOID,  NVRTDX) , TY(NV0ID,  NVRTHX) ) are  the 
coordinates  of  the  void  tracer  which  is  at  the  common  point. 

As  the  void  is  closed,  NVRTEX  is  redefined,  identifying 
another  free  surface  particle  as  the  vertex  of  the  void 
closing  region. 

8.4.2  Criteria  for  Closing  a Void 

The  void  is  closed  if  at  least  one  of  the  following 
three  criteria  is  met: 

1.  The  pressure,  PVRTEX,  of  the  cell  which  contained 
the  vertex  particle  at  the  start  of  the  cycle  is 
positive,  and  the  angle  between  the  two  surfaces  is 
less  than  .2  radians. 

2.  The  two  free  surface  interfaces  have  crossed.  (This 
occurs  only  if  the  tracer  particles  have  become  too 
sparse  in  that  region,  or  if  the  plugging  option, 
which  invokes  special  noncontinuous  rules  for  moving 
tracer  particles,  is  being  used.) 

3.  The  particle  next  to  the  vertex  on  one  free  surface 
is  less  than  1/10  of  a cell  diagonal  from  the  other 
free  surface. 

8.4.3  Method  for  Closing  the  Void 

VDCLOS  defines  four  points  which  are  uSed  in  testing  the 
criteria  listed  above  as  well  as  in  the  void  closing  procedure. 
Referring  to  Figure  8.1,  (TXl,  TYl),  (TX2,  TY2) , (TX3,  TY3) 
are  the  void  tracers  indexed  by  NVRTEX- 1,  NVRTEX  and  NVRTEX+1, 
respectively  (i.e.,  (TX2,  TY2)  = [TX(NV0ID,  NVRTEX),  TY(NV0ID, 
NVRTEX)])  . The  fourth  point  (TXP,  TYP)  is  defined  by  VDCLOS 
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such  that  the  line  through  (TX3,  TY3)  and  (TXP,  TYP)  is 
perpendicular  to  the  interface  line  through  (TXl,  TYl)  and 
(TX2,  TY2). 


If  one  of  the  three  criteria  listed  above  is  met,  the 
following  five  steps  are  taken  to  close  the  void.  These 
steps  are  described  with  reference  to  the  example  illustrated 
by  Figures  8.1  and  8.2. 

1.  The  package  B tracer  at  (TX3,  TY3)  is  moved  to 
(TXP,  TYP). 

2.  The  package  A tracer  at  (TXl,  TYl)  is  moved  to 
(TXP,  TYP)  . 

3.  The  void  tracer  at  (TXl,  TYl)  is  moved  to  (TXP, 

TYP) . 

4.  The  void  tracers  at  (TX2,  TY2)  and  (TX3,  TY3)  are 
removed  from  the  string  of  void  tracers,  and 
NMP(NVOID)  is  decremented  by  2. 

5.  Since  the  new  vertex  point  is  at  (TXP,  TYP), 
which  is  the  location  of  the  NVRTEX  - 1 tracer 
(see  step  3),  NVRTEX  is  decremented  by  1. 

Each  time  the  void  is  closed,  the  position  of  the  vertex 
point  is  changed,  and  the  void  closing  region  is  redefined. 
VDCLOS  then  applies  the  same  three  criteria  to  the  new  void 
closing  region,  closing  the  void  and  redefining  another 
void  closing  region  if  any  of  the  criteria  are  met.  This 
process  continues  until  a void  closing  region  is  defined 


* Since  PVRTEX  is  the  pressure  of  the  cell  containing  the 
vertex  point  at  the  beginning  of  the  cycle,  the  first 
criterion  is  not  applicable  if  the  new  vertex  point  is 
not  in  that  cell. 


>1  ft  i' 


Figure  8.2--  Void  closing  region  after  void  is 

closed.  The  old  positions  are  ^ 

represented  by  dashed  lines,  and 
the  original  position  of  (TXl,  TYl) 
is  represented  by  an  open  circle. 
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which  satisfies  none  of  the  applicable  void  closing  criteria. 
VDCLOS  is  then  exited  and  not  invoked  again  until  the  follow- 
ing subcycle  of  INFACIi. 


S 
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CHAPTER  IX 


ERROR  CONDITIONS 

The  following  sections  should  help  the  user  diagnose  a 
variety  of  error  conditions  which  may  occur  during  a HELP 
calculation.  The  error  messages  printed  by  tlie  code  are 
listed  alphabetically  in  Section  9.1,  along  with  the  sub- 
routines which  print  each  message  and  comments  which  suggest 
what  the  user  may  do  to  rectify  each  problem.  Section  9.2 
describes  the  optional  diagnostic  prints  the  user  may  obtain 
in  order  to  clarify  an  aspect  of  a calculation  or  an  error 
condition . 

9.1  ALPHABETIC  LIST  OF  ERROR  MESSAGES 

The  table  that  follows  is  an  alphabetic  list  of  the 
error  messages  printed  by  the  HELP  code.  Each  error  is 
classified  as  fatal  (F)  or  nonfatal  (Nf);  only  fatal  errors 
cause  execution  to  stop.  The  subroutine  line  number  of  the 
statement  whicli  prints  the  error  message  is  indicated  in 
parentheses  under  the  subroutine  name.  The  comments  direct 
the  user  to  the  most  likely  cause  of  the  error  condition; 
however,  the  user  is  encouraged  to  pursue  the  matter  further, 
as  there  are  undoubtedly  situations  which  are  not  described 
here . 

In  many  cases  the  user  can  correct  the  error  by  redefining 
an  input  parameter  or  by  coding  changes  to  cell  quantities  or 
tracer  positions  in  subroutine  INPUT.  (See  Section  7.4.)  The 
calculation  can  then  be  restarted  on  the  dump  cycle  before 
the  error  occurred  unless  the  error  occurred  in  one  of  the 
generator  routines.  However,  if  the  error  occurs  at  the 
beginning  of  a cycle,  before  SPHASE,  HPHASE  or  TPHASE  have 


changed  any  of  the  cell  quantities,  and  an  error  dump  has 
been  made  as  a result  of  a call  to  ERROR,  the  user  can 
correct  the  error  and  restart  from  the  error  dump.  (Error 
#32  is  an  example.,) 


Error  Message 


Type 


Source 


Comment 


1.  ALL  SPECIFIC  VOLUMES  Nf  CDT 

MULTIPLIED  BY  (147) 


2.  BLANK  CARD  Nf  CARDS 

(26) 

3.  CHECK  FIRST  RECORD  OF  F INPUT 

THE  DUMP  AND  FIRST  (199) 

DATA  CARD  OF  INPUT  DECK 


ERROR  CONDITION  IN  F THETAS 

THETAS  M,  MOS,  (MSLD  (78) 

(N),  N»l,  MOS) 

ERROR  CONDITION  - F THETAS 

THETAS:  THETA,  XTP,  (115) 

XBT,  YRT,  YLF 


6.  ERROR  EXIT  - SEE  F ERROR 

STATEMENT  NUMBER (9-36) 

IN  

7.  ERROR  IN  ADDTCR.  Nf  ADDTCR 

NUMBER  OF  TRACERS  (170) 

IN  PACKAGE 

EXCEEDED  NTTOTI 

8.  ERROR  IN  DETIME  CELL  F DETIME 

1=  J-  HAS  NO  DETONATION  (131) 

TIME,  K,N« 


Normal  iteration  pro- 
cedure not  sufficient 
to  conserve  volume  of 
cell . 

A blank  card  is  in  the 
set  of  cards  read  by 
CARDS. 

Indicates  user  is  read- 
ing wrong  tape  or  has 
asked  for  wrong  cycle 
or  problem  number  on 
tape . 

Flags  (MSLD)  definedby 
PTSAV  are  in  error.  Pos- 
sibly a storage  error. 

Slipline  does  not  cross 
this  cell  and  yet  PTSAV 
has  indicated  that  it 
does.  Possibly  a storage 
error. 

Refer  to  coding  near  the 
statement  number  in  the 
indicated  subroutine. 

To  activate  ADDTCR  on 
subsequent  run,  change 
NTPMX  and  reset  NADD. 
Tracer  particle  arrays 
may  need  to  be  enlarged. 

A detonation  time  has 
not  been  calculated  for 
a cell  containing 
explosive.  Check  detona- 
tion points. 


Error  Message 


Type 


Source 


Comment 


ERROR  IN  DETIME 
ROUTINE,  SEARCHING 
FOR  K-  KW  WAS  IN- 
CREASED TO  A VALUE 
.GT.  KMAX  OR  .LT.  I 


10.  ERROR  IN  INPUT  CARD,  F 
IDET= 

INPUT  CARD  IS 


11.  ERROR  IN  MFLAG  DETER-  F 
MINATION  1=  J=  K= 


12.  ERROR  IN  PTSAV  I,J,  F 
Ml,  M2,  MT,  THETA (MT) 


DETIMfi  Indicates  the  search 
(425)  for  the  K-th  cell  is 
outside  the  grid.  Can 
result  from  machine 
round  off  error,  giving 
an  erroneous  value  of 
I or  J . 

DETIME  Indicates  the  type  of 
(444)  initiation  point  (IP) 
read  is  in  error. 
Attempted  to  read 
primary  IP  data  into 
the  secondary  IP  array. 
Input  cards  are  probably 
out  of  order. 

FILGRD  Usually  a sign  that 
(134)  MFLAG(2)  was  defined 

incorrectly  in  the  in- 
put deck.  May  also  be 
an  error  in  material 
interface  definition. 

PTSAV  Slide  line  flag  (MSLD) 
(35)  for  cell  I,J  improperly 
defined.  Probably  a 
storage  error. 


13.  ERROR  IN  VOLUMES  1=  Nf  FILGRD 
J=  MFK=  VSUM=  VCELL=  (31) 


Sum  of  material  partial 
volumes  do  not  equal  the 
cell's  volume.  Can  prob- 
ably ignore  if  error  is 
very  small  and  occurs  for 
only  one  or  two  cells. 

Can  be  a result  of  inter- 
faces being  incorrectly 
defined  by  the  input  cards. 
Check  especially  that 
cards  defining  tracer 
particle  positions  for 
different  materials  adjoin- 
ing a common  interface 
segment  have  the  same  end 
points  and  define  the  same 
number  of  tracers. 


ERROR  ON  PRECEED- 
ING  DATA  CARD 


EVACUATION  OF  NIATE- 
RIAL  N I,J,N,M,ML, 

MB  TFLUX,  XMASS(N,M)  , 
SAMMP(N.M),  SAMPY(N.M) 


FRACS  DETECTED  AN  IN- 
COMPLETE SUBPACKAGE  OF 
MATERIAL  PACKAGE 
CHECK  REMAINING 
PA^  'IA.GES  Ml,  M2,  NN- , 
TXi,  TYl,  TX2,  TY2  = 

IMAX=  WHICH  IS  NOT 
AN  EVEN  NUMBER.  THE 
GRID  WAS  NOT  REZONED 
IN  THE  X-DIRECTION 

IMPROPER  GRID  DEFINI- 
TION IN  SETUP. 

IMAX  = I *= 

IMPROPER  GRID  DEFINI- 
TION IN  SETUP. 

UMAX  - J = 

INSUFFICIENT  AMOUNT 
OF  MIXED  CELL  STORAGE 


21.  UMAX  * WHICH  IS  NOT  AN 
EVEN  NUMBER.  THE  GRID 
WAS  NOT  REZONED  IN  THE 
Y-DIRECTION 

22.  MASS  EVAPORATED  DUE  TO 
ROUND-OFF  IN  TPHASE 
I,J,MFLAG,AMX,T.E. 

23.  MASS  EVAPORATED  DUE  TO 
ROUND-OFF  IN  TPHASE, 
I,J,N,MFK,XMASS,SIE,RHO 


CARDS 

C28) 


DMADJ 
(126\ 
I 145/ 


FRACS 

(562) 


REZONE 

(27) 


SETUP 

(149) 

SETUP 

(151)- 

NEWFLG 

(6) 


REZONE 

(156) 


TPHASE 

(672) 

TPHASE 

(700) 


Probably  numbers  are 
punched  in  the  wrong 
column  on  a card  read 
by  CARDS. 

Not  an  error  condition 
unless  XMASS(N,M)  is 
relatively  large  for 
evacuation  into  a neigh- 
bor cell.  Can  be  caused 
by  interfaces  crossing; 
look  for  error  31 . 

NMP(N)  not  properly  de- 
fined or  dummy  endpoint 
not  added  when  tracers 
were  generated  or 
manually  changed. 

Grid  must  have  an  even 
number  of  columns  to  be 
rezoned  in  the  X- 
direction. 

User  did  not  define 
IMAX  values  of  DX. 

User  did  not  define 
JMAX  values  of  DY. 

Code  wants  to  generate 
more  than  NMXCLS  inter- 
face cells.  Redimension 
material  arrays  and  re- 
define NMXCLS. 

Grid  must  have  an  even 
number  of  rows  to  be 
rezoned  in  the  Y- 
direction. 

Mass  of  a pure  cell  is 
too  small.  This  should 
not  occur  frequently. 

Mass  of  material  N in 
an  interface  cell  is  too 
small  or  could  not  be 
evacuated.  Not  an  error 
condition  unless  XMASS 
is  relatively  large. 


9-4 


Error  Message 


T 


Type 


Source 


Comment 


Nf 


26. 


MIXED  CELL  MASS 
TRANSPORT  ADJUSTED 
TO  PREVENT  OVER- 
EMPTYING  I,J,M,N, 
SAMPY , SAMMP , SAMMY , 
SGAMC 

NEW  SLIP  ENDPOINT 
NOT  DEFINED  L, SLOPS, 
SLOPT,TXCL,TYCL,TXl , 
TYl,  TX2,  TY2=  ERROR 
IN  PLGTCR 


NEW  SLIP  ENDPOINT  F 
OUTSIDE  GRID  SLPNDX  , 
SLPNDY, ALPHA, TXCL, 

TYCL.  ERROR  IN  PLGTCR 


DMADJ 

(59) 


PLGTCR 

(118) 


PLGTCR 

(23) 


27. 

NUMBER  OF  PARTICLES 

F 

TSETUP 

i 

HAS  BEEN  EXCEEDED 

(156) 

f 

i 28. 

NUMBER  OF  SLIDE  CELLS 

F 

PTSAV 

i 

EXCEEDS  INPUT  VALUE 
OF  NSLD 

(69) 

' 2Q 

1 

. 

REZONE  ROUTINE  CAN 
NOT  HANDLE  PLUGGING 
PROBLEMS 

F 

REZONE 

(9) 

' 30. 

SUM  OF  FRACS  NOT 

Nf 

FLGSET 

[ 

EQUAL  TO  TOTAL  AREA 
I,J,TAU,XDY2PI  = 

(191) 
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Not  an  error  condition 
unless  it  recurs  for 
the  same  cell  on  many 
successive  cycles. 


Plug  surface  extended 
into  free  surface  cell 
but  can't  find  an  in- 
tercept with  back  of 
target.  Check  tracer 
particles  of  both 
packages . 

Plug  surface  (slipline) 
can  not  be  extended  to 
a point  outside  the  grid. 
Check  the  angle  ALPHA. 

Check  input  defining 
tracer  particles  as  well 
as  value  assigned  to 
NTPMX. 

Redimension  slide  arrays 
and  redefine  NSLD. 

Plugging  calculations  do 
not,  in  general,  require 
rezoning . 

Interfaces  may  have  crossed 
due  to  tracers  becoming 
too  sparse.  If  print 
occurs  only  once  for  a 
given  cell,  the  error  is 
caused  by  roundoff  and 
can  be  ignored.  Otherwise 
back  up  and  move  tracers  i 
manually  in  INPUT  or  have 
ADDTCR  add  tracers  in  that 
region.  (See  Section  7.4.)  * 

(This  error  can  also  occur 
if  the  interfaces  are  in- 
itially defined  incorrectly.) 
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Error  Message 


Type 


Source 


Comment 


31.  TROUBLE  DEFINING 
MATERIAL  OF  CELL 
THAT  HAS  BECOME 
PURE  IN  FLGSET, 
I,J  » 


32.  TROUBLE  WITH  PRES- 
SURE ITERATION  I = 
J = ITC  = PAV  = 


FLGSET  This  can  be  a result  of 

(145)  interfaces  crossing. 

Look  for  error  31.  It 
also  can  result  from 
an  isolated  mass.  The 
user  can  "evaporate" 
this  mass  manually  in 
INPUT  on  a restart. 

CDT  This  can  be  a result  of 

/142  \ too  small  a value  of 

\229  / PMIN,  or  it  can  be 

caused  by  a discontinuous 
equation  of  state  or  one 
which  has  a negative 
slope  in  the  P-V  plane. 


I 


9.2 


INTER  DIAGNOSTIC  PRINTS 


Additional  diagnostic  prints  can  be  generated  by  the  user 
by  appropriately  defining  the  input  variable,  INTER.  The  table 
below  indicates  the  nature  of  the  diagnostic  prints  for  each 
value  Ca  prime  number)  of  INTER.  If  the  user  wants  more  than 
one  of  the  prints,  set  INTER  equal  to  the  product  of  the 
corresponding  prime  numbers.  For  example,  if  total  energy  is 
not  being  conserved,  the  user  might  set  INTER  = 5 x 7 to  see 
which  cellCs)  are  responsible  for  the  error.  In  order  to  * 
also  see  what  the  cell,  quantities  are  before  and  after  each 
of  the  phases  the  user  would  set  INTER  =3x5x7. 


VALUE  OF 
INTER 

PRINTS  FROM 
SUBROUTINE 

NATURE  OF 

INFORMATION  PRINTED 

2 

CDT 

. 

All  intermediate  steps  in  pressure 
iteration  for  each  multimaterial 
cell  in  active  grid. 

3 

EDIT 

An  EDIT  print  of  all  cell  quanti- 
ties after  both  SPHASE  and  HPHASE, 
as  well  as  the  normal  EDIT  follow- 
ing CDT  (on  EDIT  print  cycles  only). 

5 

TPHASE 

Mass  transport  terms  at  each 
boundary  of  every  cell,  plus  total 
energy  in  grid  after  each  cell  in 
active  grid  is  updated  in  TPHASE. 

7 

SPHASE 

. 

Total  energy  in  grid  after  each 
cell  in  active  grid  is  updated  in 
SPHASE. 

11 

UVMOD 

Post-TPHASE  material  velocities  in 
slip  cells  before  and  after  the 
slipline  conditions  are  invoked. 

13 

UVCALC 

Post-HPHASE  material  velocities 
in  slip  cells  before  and  after 
the  slipline  conditions  are 
invoked. 
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NOTH:  The  user  should  be  aware  that  invoking  several  of 

the  above  options  and  running  a calculation  for  many  cycles 
can  generate  a prodigious  amount  of  output,  particularly  if 
there  are  many  mixed  cells  and/or  a large  grid.  Frequently, 
especially  in  the  case  of  pressure  iteration  problems  or  an 
energy  error,  it  is  necessary  to  invoke  the  diagnostic  prints 
for  only  one  cycle  in  order  to  get  the  necessary  information. 


Vv. 


(The  reverse 
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